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ABSTRACT 


Ray methods of generating synthetic seismograms have 
been favored in seismology over the more exact wave solu- 
tions since they may be applied to problems with more com- 
plex boundary conditions. 

The ray method and code generating scheme of Hron 
(1972) as it applies to vertical incidence has been 
reviewed in detail. A new method of generating the ray 
code based on the concept of permuted partitions has been 
developed in an attempt to reduce the high computer storage 
costs associated with the generation scheme used by Hron. 
Restricting the partitions according to the number of parts 
forms the basis of attempts to further reduce costs by 
eliminating uninteresting phases. The effect of including 
geometric spreading in the amplitude computations has been 
examined and found to be of significant importance. The 


a order multiples is briefly described and 


ray method of n 
compared with the method of Hron. Finally the complete 
wave solution for a sixteen layer model is presented and 


compared with previous results. 


iv 


=o evan Hg anatdors 3 sottage ot a a te any 
“sanotsthn09 Yrs bavoc . 5 
ek 0 smadtae ontdersmae abo ‘boa boston yer aie ras is 
need 26d ssvabtoni Totty oF 2otlags at eet er | 
(ar ont ontasiensg to bodvon won A ir Thesab: at bi stven 
\Jnaad esf enoks ti4eq: retina 0 $qone3. odd fio band 9bo: 
spsrot2 iotuqmo fiptt att soubsy ot tometts ns a hing oh | 
.novk yd been omotioe norrs1909p end id tw bedstoozes ateoo 
eiv6q to 19dmun 9nd oF pnt brosos, epotstinag ons pnttotiseeh 
ud. 2t200 aoubaq ner o¢ etqmetts to etesd ant mvt 
patbulant to sastts sat -2g26dq entiesratatne pntisatmthe 
nsed zed enotrstuamoo shutitams edt wi patbsarqe atxJomesg 
edT ,sonstrogmt insotttipte to 9d oF bnvot brs bentmsxe 
bis bedtis29d ultetnd et aatats ton yob0 NJ, to bodtem 7] 
alat qmos aut yitent] 907k +o bottem edt ditw bavsqmas 4 
bs batneesyq 2 febom fsyst fsatxt2 6 yot noftuloe syvaw | 
.2dfvess 2uotverq agtw bovaEqmoed 


ACKNOWLEDGEMENTS 


Thanks are due to my supervisor Dr. E.R. Kanasewich 
for his guidance throughout the project and to Dr. F. Hron 
for his useful discussions and encouragement. 

I would also like to thank Mr. Claude Athias for the 
many late-hour discussions on topics of current interest, 
as well as my other associates and friends, Mr. Daniel Au, 
Dr. Allan Bates, Mr. Dave Ganley, Mr. Jens Havskov, Dr. Roy 
Hibbs, Mr. Joseph Lee, Mr. K. Sprenke, and Mr. Larry Sydora, 
all of whom contributed to making my stay more enjoyable. 

Special thanks are due to Mr. Dave Ganley for making 
available his computer program for calculating the complete 
wave solution to the sixteen-layer model. 

In addition I would like to thank Chevron Standard 
Limited who provided me with support in the form of a 


graduate fellowship. 


a +07 ane buat nats 
<dearsint Hove ‘Yo zafans no enoteauo ade ies nes 
euA Potne@ 7M <abrisia? bas parent un 26 trem a‘. 
woh 40 .votevel ent .0M vat hed. oven AM, 28t66 watt 1a 
evobee wives) AW bas .o%ndeqe oH “4M (sod fgseob .4M ,eddtH 
/ Jatdeye hits ovom ves2 Yn pattem of betudtitaos modw to Ths 
“Qube Tot ystrisd sved .oM of sub S15 ednsse tatosqe2 
pbstqnos ony ents shuotss 197 meyporg astugmoo 2td stdslteves 
isbow evel naeoxte ot of nobautoe avew 

 bunbaes? nowadid iiedd 8 stil | bruow T notatbbs al 
te Siete nant 0 42 dd tw oh Debiveda onw bestmed 
| .dttawol fat sxsubate 


CHAPTER I =" THEPRAY METHOD.OF HRON. (1972) 


] 
] 


ie 
lhe 
ie 


CHAPTER II 


#| 
5 


| 


Ta Bee O re CONTENTS 


Lee GC UGCA OTe mre oe Sh OREM 2G 8 ae alk ep taicars, eon as 
Ba Sac mid a CuDD. | CSeem. seats oS: Fe RR ls: 
RSL Mei enc RGM Td: 0 CIS! Meer. RE. eee ae, Urn, eine 
DAV HOME meee LIC i) CUS emerewey Liters e, Ar 0eNe craw (oe, Rem, vammnc AM at obs 


Hron's Selection Criteria and Algorithm 
nOGetnesGcencra vronyot the Ray Code... «24, 


- THE PRACTICAL APPLICATION OF THE .RAY 
PUNT MOBS kOe wT tet) EL (AS 7k) Rs Any eed canta ane al ar ie! 


Tes Mat rg xa ee elec eie eee cae en oo ee Bae Seca ee 
Partitions, Compositions and a New 
Algorithm for the Generation of the 
Kae Me lier C OLUEm 7 carn lematareralar ee Tarra. (arene mya vss ee, 
Lies aay Coto NIROms COU Ciiiers © emia caer Polen ads. 


Lie Bermnutat Onis Ub FOU DI Gitta ore 5s eee clea 


Cost Comparison Between Hron's Algorithm 
and= the Permuted= Pant tion) Approach an... s. 


Reg ra Crue Cama Tele CO TIMLM Cte catia Sine tan arenes Se 


Synthetic Seismograms with Geometric 
SDE Ch NG arcetee ot ay cates tame a aa he ner era a coeie gn vette Mag once aes 


CHAPTER @ TUL, HESRAYS MET ODe OF ntPorDER 
Pt (aT OL PL Osten eer be dtc aati: nr a nase (omnss ler wear 
Bia’ SCOP GO ILC OWES cancer tieetsacdy. so aly a cere onan ope te yer cy aes 


o. 
oe 


] 
2 


The Coding Scheme and Method of 
GeTIeWatl ONL teenies oer ec tarice ein. c heen tere St Mette whanel cscs 


vi 


ne 


errs tr iro Ci te perme Pray 
ee See oe eae eh ers 2poteih Speman ht B« f : ’ 
rie eee Woe eas ie. een ohasag: ; ee be 


Pe aS er pom. alge a sant ot Ah 
Mee). Aw. Riker hog.’ restyann an spare 2a) ~ LE a TSA 
fis 


Cite daist oth FS) 


ee G@is'es © be twee 78 * 


| “pm Seek mar eeth $.$ 


aS cet 
BS eae ON Ad | a aniiiiibaeltwe oft €.9' 
a3 casas sting ~aipuatine nofinggatet oT FS 
uth? TA oy a mchrtqiod te09 8.8 * | 
gf ee ee. he ee bagumret vids bos _ 
a een cesteeeseeRntaokshta6? Betortsesh 8.5 
ate 5 chinese ie attads nye £3 


ee ee tne. to C eG 9 tbh @ vhs prtanorge 


naasatn x0 oontan (ag set ~ hr SOND 


Pie 
Ja iia? bin Aa ee aes wedia een as 
r 


We te ww Maya wi ee ad eae ale .7oqernes dtescé f.£ 


eC ls: 4 nGA2M a8) Shatin? pated a4T $.8 


r, pits fee Paws Je Ee nig e+ ides € bts , fol tareped 


wv 


ie a ‘i 


Page 


3.3 Comparison of the Method of Nth Order 
Multiples and the Ray Method of Hron 


CO Or7: a) ree oe mere sr See taVer at Shotts cite he uke 03 Glare he 49 
32 4 oui CO OG EGEEMU CT DALCISe severe d.chehou be etencas eebene 54 

3.5 The Complete Wave Solution for the 
lbetayer Model at en ALeantithm fan the. 2... oh) 
SB ODER CO MICH US:O.N) Stave emer tent einer ate, cdot laltattstene es etree es atts 98 
BIBLIOGRAPHY@%7.© .A% Stee Eee ott hs labels teeter ame fo stiles oteltte oPetemieteleny« 61 
AP PENDIX Ato omtceratcle te ch ststece, stelete tcl ce sretetecs yet sc: Lela eRe a peu age es 62 
PEND UKM Ge iaonbevscc tes lalate: ohecstereysie ta \ctale voters erie etic nits ve uec lene tenet ens 64 
APPENDIX C. SHORE AAT, 5 PRR: PRMD at, OR CITR is ne ese Hate Pas Fe te feito to He Ve 65 


vii 


ao BS emt e b place) wale. 


a 
add yo? nol tuto? ave! 
bad se wie 6 b 6 > het Om ip_g eee bow 878 ya Oe Bdaed , 
- : Lice al ae = ‘ 
sul tulateh san waa den tases ek» eee a.€ 


7 . 7 


rs ONDE pan Hd eden s ReduTNy Ae nied einen sls 2a MOOR 


| 
seagate « bast rp ceshtad aoe ohad-ashies ae ae toe XTOWIT9IA 
cine bead ce W MAE Alen ell Gea ey oe - 


ee oer Peers KLOWI949A 


ens @ee @ © © 82 e@ eee FP tee? 


TABLE 


I 


14 


5 


Lome Ora PABILES 


Aba Cone OUCDUt. TOY WHMS. SLAY. =S.At ne ctee s o> > 


Pascal's Triangle - ae BEM <a Gen Weds vate tates Ge ceknire ts) 


Demonstration of an Algorithm for the 


Creation of Symbols for Kinematic Analogs.. 


Comparison of the Total Number of Kine- 
matic Analogs Mt(LAY,HMS) and the Reduced 


ROW DiMenSsOneDIMetOrehMS am EAY eet, ole cs 
Hota SpacesReguared: fOr Nill) gem wes e cis 
The Unrestricted Partitions and Composi- 
CroOns* OF GtheaNUMDENSH IN COR 4. yt rene is cyeteteteatoi 
Raped OnseOt NebDyONUMDer On RantSoe ss . cm. 
Initialization of Matrix N(k,L) and 

MECL OVE) oii te eerie etch ten nenetg Saneeetadets (onesie 
Generation of Sequence (02) o.4 .. sc 66 cigas 5 
Generacroneor SeQUEN CEN (3 nie eni 2 cts dss tare. 


Comparison of Vector Space Requirements 


for the Generation of the Kinematic Code... 


Cost Comparison of Two Different Programs 
for the Construction of Synthetic Seismo- 


CAMS m.. SMe eeNs PORN Rl ie el avousieaic om retacs feet cue 


Cost Comparison of Four Levels of 
Restriction in the Generation of the 


h6=Fayer S¥YNthetic<Selsmogram cs, ws sie toc Wnts 


Average Relative Percentage Differences 


fOK the .b6=LavereMo deuce, Mine iter neers wis vere’ 


The Average Relative Percentage 
Difference With and Without Geometric 


SD GG AGU Ciece erat cant teen hak ois nea Wee cral teu coeee war waere 


viii 


1g 


20 


21 


ee 
Za 


30 
30 
si 


32 


38 


39 


4] 


44 


, bpd eaten cists 


Z oe 
din ieee renee 0 oa = a ‘a ae Tsotgut 
ont Ae ee 0 - Site 2! faoe89 


on ad +O} mad Fn TA 16 sank woe 
.. 2pofeak et npcat eto oe Yo notise) 


-sht¥ to t9dmult fstoT sft to qoztrsqmod 
Fai eh afd bos (emi, VAL), i. 3 Tsah ottem 


a et ee Y¥AL = 2MH vot stanomba wos 
een 2 4, Ln Tot beytupe? sosqe [stot 
-~feoqmo) bas enotitirs baxaratesynu sat 
ha casa nesesd OF T eyedmun ofd to enots 
ee ig .23769 To yedmul Yd A Yo enordt31e9 
bas (4,4) xinteM to norss Phatstol 

Dalgie spielen Po Valen a ee sig rd Tt b totosV 
stihe: _pecsecevs(S) goneuge? to notoerened 
i See eka eee wea (EB) SonSUpse to nottavensa 


2inomevhupsa edsge vos99V Yo soztisqmo? 
,.9b09 attemaniX oft to noffevsnad edd 10F — 


a Se dnoyetttd awT to doztrsqmod t202 


-om2te2 attaftaye to notsauitenod edt vot 
“ee 8 @ te eweee *-? ee re Ler ee eG 


v8 ayot TO “foziveqgmo)d F209 
sna Thiet Ane ee noTsar+sear 
wll ee ogy fe 15%6 1-31 


ganna ght 6 aber Yeeeah sae nen 
Sous Sgee th 


it ahh hana a s97ge 


Ot 


TABLE 


16 


ihe 


Page 
Average Relative Percentage Differences 
Uist ouN CleOnde ce MU ti piles eee vane tre 52 
Cost Comparison of Generating Primaries 
Plus Up to Second Order Multiples for the 
[6akayerneModek. oF. Gay. Path. Seqmegns. case as 54 


ix 


ny oe 
is, ern cad mere i 


a 


hte gf oe R 
CPR Ale vs 
} Lh jerr~'e 
Thor ’ 
. 


7 
: 


FIGURE 


Ls 


14 


i) 
16 
17 


lou Ore EIGURES 


Illustration of the Concepts of Kine- 
maticpand,DynamiceAnalogSa¢as. fav.Bee.....- 


Classification of Ray Path Segments........ 
Gontinuity of RaygRaths. har cgecenee. Fae. ..- 


Distribution of Trivial Elements into 
"Empty" Elements of the First Class........ 


The Three Ways of Selecting One Trivial 
Element from Three Elements of the First 
CHAS |S tot ecedeh otencishe acta neat e cea eka eee ee Ee Tee 
The Dynamic Analogs with Code (3,3;1)...... 


Procedure for the Generation of Synthetic 
SELSMOGKANS Lees. sent. eee eee eee. 


LhesPar titi oning? Algo rachincereecs tee eet re 
Lhe Permuteatd ona ligio rat nitseenws eee ceeteieton ss. 
ReEStMictedMParn titi OMS: ae ee ote eR he 
Sixteen Laver Mode laveiwir vtec. coeee ee Caer. 


Spike Synthetics Generated Using 
Different. Levels of Res CY 1CULON = werencst.cterinre 


The Sonigram and Synthetic Seismograms 
forthe W6=Uayver ModetL Tah screen ceee eee 


A Comparison of Synthetic Seismograms 
With and Without Geometric Spreading....... 


First and Second Order Multiples..........6.- 
Two Sequential First Order Multiples....... 
Spike Sequences for Severity Level 3 and 


Primaries Plus Up to Second Order 
Mol Titd DY CS sc oc hc-coevas east oan teeias 5 arvepece rear austeier er eann 


38 


43 


45 
46 
49 


50 


a a uh Yoh ; 
My Pree  a: in wae 2 


| See tea 


7 


ee »(f28,£), sbod it tw 2pof snk eet adT 
as 2ttendtinye to nottsxened dt 10t ov 14 
Bek Vi bwisicateweaietey |. sauts P90p 2 pitnots ttv69 ait 
Bas V6sadcocwa fie wes , ‘mds tvop tA noktssumys9 ent 
Be) OY anc a ala oie a eA oN ..@norsts159 badotrteaa 
, Serre eee: tes nyefQboM 18Y81 mgsdxt2 
pnt2u betes ane Seren oat 

ig hain Sout Mofsotysz: to elaved : eh 
emsypomared attadtaye ba giT 

e*#evee ees a oF 


FIGURE 


18 


19 


20 


2] 


Synthetic Seismograms for Severity Level 
3 and Primaries Plus Up to Second Order 


MUG) OSS 5 eae s Lites baccetat 5 ee Pe aa Pe 


Comparison of Spike Sequences for Nth 
Order Multiples and Nth Order Surface 


MU CTP Vas s bY Use ta we ee beats TRE eee OE CO 


Comparison of Synthetic Seismograms for 
Nth Order Multiples and Nth Order Surface 


MG Tp | OS ine a da Tee eka eon ee awaeees UU a TAMU 


Comparison of the Unrestricted and 
Complete Wave Solutions to the 16-Layer 


MOOG) Fee PAT; I PSSRION, OD; ADS, SUNT POR 6 i 


x7 


Page 


Jd 


56 


a7 


59 


rey ¢ vt! iy Uae oF + on 
te ; pe 2 esa a 
i’ 


we sin 


viet Hagens has 
* ‘ 


Xe ee 


| ; A if A 
' j i » ric j ' 
i ies ) 
. 
Ng 
vu) >| 
A 
aay Ol 
7 ra a : 
‘ vit 
Le, ) : 
- AL og Pe on 


AAS yh os en mm 
id 


CHAPTER I 


THE RAY METHOD OF HRON (1972) 


ete tnt OUueULOon 

A common approach taken in the construction of syn- 
thetic seismograms is the ray tracing technique. In this 
method a group of rays is selected and the amplitude contri- 
bution and time of arrival of each ray within the group is 
determined at a particular location on the surface. The 
problem lies in the selection of a group of rays from the 
infinite number of rays travelling from the source to the 
receiver. There have been various approaches taken to this 
problem and these will be discussed and comparisons made 
using synthetic seismograms at vertical incidence. 

Practically, synthetic seismograms at vertical incid- 
ence are of limited interest however their construction can 
serve a useful purpose. In exploration work, for example, 
when the depth of energy penetration is large compared with 
the shot-geophone spacing, a synthetic seismogram at verti- 
cal incidence may approximate to a high degree of accuracy 
one constructed taking into account the shot-geophone spac- 
ing. This may be especially true if the near surface struc- 
ture has little in the way of velocity contrasts. Of 


greater importance though is the fact that synthetics at 
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vertical incidence offer a reasonably inexpensive method 
of comparing different ray selection criteria. The results 
of the comparison can then be used as a basis for ray path 
selection when proceeding to non-vertical modeling. 

One approach to the problem of ray selection is due 
to Hron (1972) with the aid of prior work done in the Soviet 
Union. The remainder of this chapter consists of a review 


of Hron's work as it applies to vertical incidence. 


1.2 Basic Principles 


Rays which travel from the source to receiver along 
different ray paths with identical travel times are defined 
as kinematic analogs (Hron, 1972). Within each group of 
kinematic analogs there may exist subgroups of rays which 
have equal amplitudes and which are defined as dynamic 
analogs. The necessary and sufficient condition for kine- 
matic equivalence of two different ray paths or phases is 
that each must have equal number of segments in each layer. 
It will be assumed that both the source and receiver are 
located on the surface, hence the number of segments in each 
layer must be even. That is, for every down going ray in 
layer i there must exist a corresponding up going wave. 
Failure of this condition results in an incomplete ray path. 

To describe and identify groups of kinematic analogs, 
J integers are used as the parameters, letting the symbol 


of the group be 
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where ny is an integer equal to half the segments in layer 
j and J is the total number of layers involved. 

Figure 1 illustrates the concepts of kinematic and 
dynamic analogs. The three phases drawn comprise the group 
of kinematic analogs with symbols (2,2), each ray path 
drawn having two half-segments in the first layer and two 
half-segments in the second layer. Within this group of 
kinematic analogs there exists a subgroup composed of phases 
(b) and (c) which have the same number and type of inter- 
actions with the interfaces and which are called dynamic 
analogs. The modification of the coding in (1) to include 
a description of the dynamic properties of a particular 


phase will be dealt with later. 


Figure 1. Illustration of the Concepts of 
Kinematic and Dynamic Analogs. 


(a) (b) (c) 
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For a more complete discussion of kinematic and 
dynamic analogs several important concepts are necessary. 
A coupled segment is defined as one down going and the 
nearest up going segment in a particular layer. Since the 
number of segments in each layer is even there are ns 
coupled segments in each layer j. For example in Figure 2 
segments 2 and 3 constitute the first coupled segment in 
layer 2. Segments 4 and 7 constitute the second coupled 
segment in this layer. 


th 


An element of the j class is defined to be the con- 


tinuous chain of segments bounded by a coupled segment in 


the jth 


layer. Any ray path bounded by a coupled segments 

in layer 2 for example will be defined as an element of the 
second class. Hence segments 2 and 3 constitute the first 

element of the second class and segments 4, 5, 6 and 7 


constitute the second element of the second class. Since 


Figure 2. Classification of Ray Path Segments 
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th 


each element of the j class is bounded by a coupled seg- 


ment they cannot overlap and since there are ny coupled 


segments in the aa layer there are exactly n, elements of 


th class 


the gth class for all j. If an element of the j 
consists of only two segments such as the first element of 
the second class consisting of segments 2 and 3 it is called 


th class is 


a trivial element. A normal element of the j 
defined as consisting of at least one coupled segment in 
the (j+1l)st layer. An example would be the second element 


of the second class consisting of segments 4, 5, 6 and 7. 


1.3 Kinematic Analogs 

Once the kinematic code is specified thus establish- 
ing a particular group of analogs the segments must be 
connected to form a continuous ray path. Hron (1972) 
defines continuity in terms of elements of the eon class. 
He states: 

"The required continuity of the chain of segments 

is satisfied if at least one element of the jth 

class j = 1, J - 1, is normal so that the segments 

in the jth and (j+1)th layers are connected. More- 

over, all elements of the (j+1)th class within the 

element of the jth class must be linked together." 

These criteria fail to ensure continuity as can be 
demonstrated by referring to Figure 3. 

In Figure 3 there is one element of the first class 
consisting of segments 3 to 10. The two elements of the 


second class involve segments 1 and 2, and segments 4 to 9. 


Finally there are two trivial elements of the third class 
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Figure 3. Continuity of Ray Paths 


consisting of segments 5 and 6 and segments 7 and 8. The 
element of the first class and the second element of the 
second class are both normal. Moreover, the element of the 
second class within the element of the first class and both 
elements of the third class within the element of the 
second class are linked together. Despite satisfying 
Hron's criteria the ray path is not continuous. 

In order to ensure continuity at least one element 


of the th 


class, j = 1, J - 1 must be normal. In addition, 
the (K+1)st coupled segment (if it exists) in layer j must 
be traceable from the start of the ya coupled segment in 
the same layer. 

To determine the total number of different ray paths, 
Ni (my sMoos--oMy) which may be constructed given an seg- 
ments in each layer j = 1,...,U0, one must always start from 


ny elements of the first class which must be linked together. 


In the case of only one layer (J=1) there is only one way 
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of connecting the nj trivial elements of the first class 

so that N,(n,) = 1. If the number of layers is two (J=2) 
with 2n, segments in the first layer and én, segments in 
the second layer, the total number of different kinematic 
analogs is given by the number of different ways of distri- 
buting the no trivial elements of the second class among 
the ny elements of the first class. An example may serve 
to illustrate. There are three ways of distributing the 
two trivial elements of the second class among the two ele- 
ments of the first class drawn in Figure 4. These three 
different possible associations are shown in Figure l(a), 


(b) and (c). In Figure 1(a) one trivial element has been 


Figure 4. Distribution of Trivial Elements into 
"Empty" Elements of the First Class. 


distributed into the first element of the first class; the 
second trivial element going into the second element. In 
1(b) both trivial elements have been distributed into the 


second element, the first element being left "empty". In 
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l(c) both trivial elements have gone into the first element 
and in this case the second element has been left "empty". 
As can be seen the problem is identical to the distribution 
of balls into pockets, some of which may be left empty. 
According to Appendix A the total number of different ray 
paths possible given én, segments in the first layer and 


2n. segments in the second layer is 


For the three-layer case the number of different ways 
of distributing the N3 trivial elements of the third class 


into each arrangement of the no elements of the second 


Class 1S 
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But there are C, different distributions of the No 
2 


elements of the second class among the ny elements of the 
n,tn,-1 notn,-1 

first class. Hence there are C, . C. 

2 3 

possible kinematic analogs given ny elements of the first 


different 


class, uP elements of the second class and nN trivial ele- 
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ments of the third class. 
In general, the total number of different ray paths 
belonging to the group of kinematic analogs (Ny sMos-+-oMy), 


J > 2 is given by 


1.4 Dynamic Analogs 


As defined earlier, dynamic analogs are groups of 
rays that arrive at the receiver at identical times and 
with identical amplitudes. As was shown in Figure 1, a 
particular group of dynamic analogs exists as a subset of 
a specific group of kinematic analogs. Therefore the kine- 
matic code, (NyoNos-++ ony) must be identical for all ele- 
ments of any particular group of dynamic analogs. In addi- 
tion to sharing the same kinematic code, each element must 
be reflected mi. times, j = 1,...,J0, from the pen interface. 
This ensures equality of amplitude. 

For a complete description of a particular group of 
dynamic analogs (2J - 1) integers will be used. The first 
J integers will consist of the kinematic code which along 
with the velocity structure of the layered media determines 
the time of arrival at the surface. The kinematic code is 
followed by J-1 integers (my sMy5-++sm,_7) which represent 


the number of reflections from the jth interface. The com- 
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plete code as shown in (5) allows the amplitude of the 


phase to be determined. 


(Ms Nos wees My5 mys My, «.-s My_y) (59 


It will be noticed in (5) that the number of reflec- 


GrvOnS. My> from the gth or last interface, has not been 
specified. The reason is that my = Ny; To understand this, 
th 


note that no normal element can exist in the J layer. 


The elements therefore, must be trivial. It can also be 
stated that for continuity, a trivial element must consti- 


tute a reflection. Since there are nj half-segments in the 
gth layer, there are ny coupled segments which form nj tri- 
vial elements and ny reflections. Hence ny) = ™)- 

The values which each m. can take on is dependent on 
the values of oe and Nea? j = 1, Jd - 1. To investigate 


the nature of this dependence recall that each reflection 


th 


in layer j demands a trivial element of the j class. 


h 


Recall also that at least one element of the zc class must 


be normal. Hence for ns elements of the jth class, at most 
nj - 1 can be trivial and there can be at most ns - |] 
th 


reflections from the j interface. That is 
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Now consider two cases: 


Case 1: ea Pies 5 


th class are normal then 


th 


If all n, elements of the j 


the number of reflections ms from the j interface must 


equal zero. However it is only necessary that one element 


th 


of the j class be normal so that for this case 


3 
[Vv 
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Case 2: on <n 


ipa j 


th 


In this case not all n, elements of the j class can 


be normal. The greatest allowable number of normal elements 


of the righ class is Nas: Hence there must be at least 
(n, ~ N41) trivial elements of the jth class. For-this 
case then 

m; > (n, - M541) (8) 


Collecting (6), (7) and (8) yields 
= clas Metal fs (9) 


MAX(0, n, - N47) <m 


Using (9) it is possible to generate the symbols for all 


dynamic analogs given the kinematic code. 
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The number of dynamic analogs may be defined as the 
number of ray paths which satisfy the conditions inherent 
in the ray code, (n,,n5,...,M)3 my),m5,...,m)_,). In order 


h 


to calculate this number consider first only the ie and 


h 


(J-1)st layer. The n, trivial elements of the iy class 


must be distributed among the n elements of the (J-1)st 


J-1 
class leaving M3 @] of them trivial. The solution is given 
by the product of the number of different ways of choosing 
my_; trivial elements from the n)_, elements of the (J-1)st 
layer, and the number of ways of distributing the my trivial 


elements of the gth 


class among exactly (nj) - m1) normal 
elements of the (J-1)st class. The problem is equivalent 

to the distribution of m balls into (n5_j - my _4) pockets, 
none of which may be Sipe, The result which is derived 


in Appendix B is given by 


(10) 


If the same reasoning is applied to the remaining 


layers the expression for the total number of dynamic ana- 


logs with the ray code (nj,N55...5N)95 My sMyoe++2My_1) iS 


asby0 nT Liane piles sek 
ine ch intes pheno! Ned 
geetn “Sy: itt to zanomata istre yn eer | 
#e(T-b) ody Yo, atnsnais yop ad onons baduarazeth: od teum 
vovio. at nottufoe sAT  .fetytns sit Sep! rat tnlvsot eesto 
entzoots to 2yew basi tb vo vodmun oft Yo soubovg oft ud 
te(f-b) sdt Yo etnomals [- ay ons movt atnamete tetvtasd 7 y™ 
Patviad pm edt onttudtate kb 40 246M to voi ill bre .voysl 
femvon Cr ops pM) yitosxa pivanis ‘2eBTD %% site to etnsmels 
tnahevhups: et mofdovq oat .ee6ts a(t =0) ont to ednamate 
2491909 (pop - plgm) OFnt et tsa uit 0 ndtauaing atb oft of 
baviveb 2t dotnw $tw2e% ent. tons a vse oid ‘to snon 
hy ehh rh nats er 8 nites at 


L(t + pit) 


RTT ae he Sahl 4 


“TN 


. (at) Snes, A ne x 
| Pape ett 


' a ' : } F 


iinsals intiile et pnénoenst Givi ond 41 


: a F ee 7 
a am ik a mn 
Roa pean oe teal? eee, 
¢ - m 
; . . 


| ri ros bart tenor ous 108 notaserqKs: afd diese 


the) 


(11) 


A simple two-layer example with ray code (3,3;1) may 
serve to illustrate the procedure. From the ray code it is 
observed that there are three elements of the first class, 
one of which is trivial and three trivial elements of the 
second class. There are three different ways of choosing 
the single trivial element from the total of three elements 


of the first class. These are shown in Figure 5. 


Figure 5. The Three Ways of Selecting One 
Trivial Element from Three Ele- 
ments of the First Class. 
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The three trivial elements of the second class must now 
be distributed among the normal elements of the first class. 
For each of the three cases in Figure 5, there are two 


: 3 
possible ways of doing this for a total of C5 Cy = 3-2 = 6 
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dynamic analogs. These are shown in Figure 6. 


Figure 6. The Dynamic Analogs with Code (3,331) 
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The advantages of the classification and coding scheme 
used by Hron should now be obvious when dealing with plane 
layers. Instead of calculating the amplitude contribution 
and time of arrival of each of the phases shown in Figure 6 
separately, the amplitude and arrival time are calculated 
only once. The amplitude value when multiplied by the 
number of dynamic analogs as determined from (11) yields 
the total dynamic effect of the group. In terms of the 
number of computations, the savings grow as_ the number of 
segments increase to accommodate a greater number of layers. 
In such cases the number of dynamic analogs for a particular 


ray code can easily reach several hundreds. 


1.5 Hron's Selection Criteria and Algorithm for the 
Generation of Ray Codes 


Before designing an algorithm to generate ray codes 
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as described in the previous sections it is necessary to 
select a family of rays from the infinite number of possible 
choices. Hron chooses that family of ray paths such that 

a maximum of HMS half-segments are allowed in any one phase 
Where HMS is a specified input parameter. Table 1 shows 
typical output using this criteria for HMS = LAY = 4. LAY 


represents the number of layers in the system. 


Table 1. Typical Output for HMS = LAY = 4 


Group Number Ray Code 

1 (Hp) 

2 ae 

3 eles 0 

4 (3) 

5 (sees HSMP) 

6 Clix 3.0) 

7 ( 1 sAll 3.19510, $10)) 
8 (4) 

9 Coys) 
10 (252.30) 
11 (2 1) 

LZ Ci Gay 

ics! (205%, 0) 
14 aaa 
15 P2050) 
16 Oa RAG) 


If HSG represents half the number of segments in a 
particular group of ray paths, then HSG can take on values 
from 1 to HMS. These HSG half-segments can be distributed 
among the top LM layers where LM = MIN(HSG,LAY). The restric- 


tion on LM is to ensure continuity. If LBT is the number of 
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the deepest layer penetrated by the ray with HSG half seg- 
ments, there are as many different groups of kinematic 
analogs as there are possible distributions of HSG coupled 
segments among LBT layers leaving none of them empty. 


According to Appendix B this number is given by 


M(LBT,HSG) = HSG - 1)! =a 


HSG-1 
BT > T)I(HSG - LBT)1 (12) 


LBT-1] 


In accordance with the selection criteria, LBT is allowed 
to vary between 1 and LAY and HSG varies between 1 and HMS. 
Therefore the total number of different groups of kinematic 


analogs is given by 


HMS LM 
M,(LAY,HMS) = I y  M(LBT,HSG) 
HSG=1 LBT=1 
per pt kfenmet is 13) 
HSG=1 LBT=1 % 


. .HMS=1 _ ,HMS-1 
In (13) consider the term Ch M_q = CMIN(HMS,LAY)-1 


where LBT = LM and HSG = HMS. Since the condition 
HMS - 1 > MIN(HMS,LAY) - 1 must hold in order for the term 
to have a non-zero value it may be rewritten as tins 


where LMX = MIN(HMS,LAY). Expanding (13) now yields 
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See ha ab GOs EA ak 
we in oe tinee Ea ee 
LMX-1 . .LMX-1 , ,LMX-1 
* Crux-3 * Somx-2 * SoMx-1 (14) 


Recalling that Pascal's Triangle tabulates the values 
of Cl, it may be noted that if n = HSG - 1 and r = LMX - 1 
then (14) represents the summation of up to LMX values in 
each of the first HSG rows of Pascal's Triangle. For 
example in Table 2 the sum of the values within the boxed 
area represent the total number of kinematic analogs in six 
layers with the maximum number of half-segments, HSG, being 
equal to eight. 

In order to generate the codes for the ray paths 
governed by the selection criteria, a matrix N(I,L) is set 
up, each row containing a different kinematic code. The 
column dimension must be set equal to MIN(HMS,LAY) but 
although the total number of kinematic analogs is given by 
(13) it is not necessary that the row dimension assume this 
value. Instead, a reduced row dimension may be substituted. 
Consider Table 3 which was obtained from the last eight 
rows of Table 1. The last column in Table 3 with the head- 
ing M(LBT,4) is seen to correspond to the fourth row of 
Pascal's Triangle. The three rows in Group II were all 
obtained from row 8 by subtracting integer values 1, 2 and 


3 from the first column and adding them to the second column. 
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Similarly row 12 was created from row 10 and row 13 from 
row 11 by subtracting 1 from the second column and adding 
it to the third. Row 14 was created from row 11 by sub- 
tracting the value 2 from column two and adding it to 
column three. Finally row 15 was obtained from row 14 by 
Subtracting 1 from column three and adding it to column 


four. 


Table 3, Based on Hron (1972). Demonstration of 
an Algorithm for the Creation of Symbols 
for Kinematic Analogs. 


Row Number Row Vector LBT M(LBT,4) 


8 & - I ] ] 
9 3,1 

10 Ene - Il 2 3 
11 1,3 

12 2,1,1 

13 1,2,1 - II! 3 3 
14 Tegelegee 

15 1T,1,l,] - IV 4 1 


Since group (i) may be generated from group (i-1), 
i= II, IV, in the example shown, the maximum number of 
rows necessary in the matrix N(I,L) is equal to twice the 
maximum value of M(LBT,4). In this case the value equals 


2* 3 = 6. In general the reduced row dimension DIM is 


given by 
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DIM = 2 * MAX(7,—\"MS eC ) (15) 


where 
Lis) lL. < MINCHMS:, LAY): 


A comparison between the actual number of kinematic analogs 


and the reduced row dimension DIM is given in Table 4. 


Table 4, Comparison of the Total Number of 
Kinematic Analogs Mt(LAY,HMS) and 
the Reduced Row Dimension DIM for 
HMS = LAY. 
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CHAPTER II 


THE PRACTICAL APPLICATION OF THE RAY 
METHOD OF HRON (1972) 


2.1 The Matrix N(I,L) 


In the last chapter it was demonstrated that by using 
the concept of dynamic analogs a saving in time could be 
realized when calculating the amplitude contributions of 
the various phases. Despite this advantage there is a 
serious limitation in the algorithm for the generation of 
the kinematic code as presented in 1.5. The limitation lies 
in the amount of additional space necessary for the matrix 


N(I,L) as the number of layers are increased. Table 5 shows 


Table 5. Total Space Required for N(I,L) 


L(LAY) 1 (DIM) I-L 
2 2 : 4 
4 6 24 
6 20 120 
8 70 560 
10 252 2,520 
12 924 11,088 
14 3,432 48,048 
16 12,870 205,920 
18 48,620 875,160 
20 184,756 3,695,120 
29 705,432 15,519,504 
24 2,704,156 64,899,744 
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the total amount of space (I-L) needed for the matrix 
N(I,L) when HMS = LAY, LAY = 2(2)24. From the table it 
can be seen that the amount of space required for the 
matrix N(I,L) grows very quickly. It is clear that if a 
model consisting of many layers is to be used, a different 
approach would be useful to generate the kinematic code 


more efficiently. 


2.2 Partitions, Compositions and a New Algorithm for 

the Generation of the Kinematic Code 

A partition is defined as a collection of integers 
with a given sum without regard to order; the corresponding 
- ordered collections are called compositions. The integers 
which form a partition or composition are called its parts 
and the number which is the sum of these parts is known as 
the partitioned or composed number. By convention repeated 
parts are abbreviated by the use of exponents. For example 
2éce 15 WhieLen as le Table 6 shows the unrestricted parti- 


tions and compositions of the numbers 1 to 4. 


Table 6. The Unrestricted Partitions and 
Compositions of the Numbers 1 to 4. 
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Restrictions can be imposed on partitions and composi- 
tions as to the kind (even or odd), size, number of repeti- 
tions, or number of parts. Table 7 shows all partitions of 


n for n = 1(1)8 grouped by the number of parts. 


Table 7. Partitions of n by Number of Parts 
from Riordan (1958). 


Number of Parts 
n Va 2? 3 4 5 6 fh ot 
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Recall from the previous chapter that the selection 
criteria of Hron (1972) entails the distribution of up to 
HMS segments among a maximum of MIN(HSG,LAY) layers. For 
each value of HSG the resulting distribution yields a 
collection of kinematic codes. The integers in each code 
from this collection must sum to HSG. Hence the collection 
of kinematic codes represents a composition of HSG with the 


restriction that the number of parts of the composition, 
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m < MIN(HSG,LAY). Equivalently the collection of codes 

may be viewed as the result of first partitioning HSG 
followed by generating the unique permutations of each 
partition. For example the composition of the integer 4 

in Table 6 is identical to the kinematic code (row vectors) 
in Table 3 of the previous chapter. Thus the basis for a 
new algorithm to generate the kinematic code is established. 
A general procedure for the generation of synthetic seismo- 


grams using this approach is shown in Figure 7. 


2.3 The Partitioning Algorithm 


In the computer program written to generate synthetic 
seismograms by this meltiod, the partitioning of the integers 
1 to HMS was carried out using an algorithm presented by 
Lehmer (1964). The flowchart for this subroutine is shown 
in Figure 8. This procedure generates those partitions 
whose number of parts, m, satisfy p < m <q. In order to 
generate all partitions of a given integer, p is set equal 
to 1 and q is assigned the value n where n is the integer 
being partitioned. In this manner all partitions of n are 
generated beginning with n and ending with 1,1,...,1. The 
vector "a" of length q stores the partition currently being 
generated and a variable T in conjunction with m ensures 
that this partition will be unique. The generation of the 
current partition or kinematic code does not require the 
manipulation of previously generated codes. This results 


in great savings of storage space as will be seen. 
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Figure 8. The Partitioning Algorithm. 
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2.4 The Permutation Subroutine 

As seen in Figure 7, after each partition is obtained 
the distinct permutations of the partition are generated. 
Each of these permutations including the original partition 
represent a particular kinematic code. There are several 
algorithms available for the efficient generation of permuta- 
tions, e.g., Paige et al. (1960), Johnson (1963), Wells 
(1961). They are all based on the fact that there are n! 
permutations of n distinct marks and also n! (n-1)-digit 
numbers of the forn: a,l! + a2! ee eee lia ue 
OPK a. < i, i = 1(1)(n-1). The permutations are generated 
by setting up a one-to-one correspondence between the n! 
permutations and the n! signatures (a)2a52---2a,_4)- Despite 
being efficient both in terms of processing time and stor- 
age costs, a totally different approach was necessary. The 
reason is that these methods are designed to produce all 
the permutations of n distinct marks. If the marks are not 
all distinct as is the case generally when dealing with 
partitions, repeated codes will be generated when using 
these methods. In order to prevent this from happening and 
still retain this approach, large amounts of storage space 
would appear to be necessary. Hence a different method of 
producing the distinct permutations of a given partition 
was developed and the flowchart is presented in Figure 9. 

The permutation subroutine makes use of a matrix 
N(K,L), the column and row dimensions both being equal to 


the number of layers (LAY) in the model render considera- 
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tion. Before beginning, the current partition is inserted 
into the last MAC positions of each row in the matrix where 
MAC represents the number of parts in the partition. The 
algorithm is constructed in such a way that all distinct 
permutations involving only the last k integers of the 
partition are generated before permutations involving only 
the last k + 1 integers are attempted, k = 2(1)(MAC-1). 
This is done through the use of a vector J(I) whose value 
indicates the column position of the last component in the 
partition that the integer in column I of row I has attempted 
to exchange places with. The value of J(I) is always 
greater than I and the interchange of integers is always 
between positions I and J(I1) in row I of the matrix N(k,L). 
In order to ensure that the permutations will be distinct 
two criteria are used. The first which should be obvious, 
is that no two integers that have the same value will be 
interchanged. The second criterion is that if J(I) > I + 1 
and N(I,J(I)) = N(I,k), k = (I+1), (J(1I)-1) then no inter- 
change is made. A simple example may serve to illustrate 
the general procedure for the generation of permutations 
using this method. 

Consider the partition 2133 which has been written 
into the last four positions of the last four rows of the 
matrix N(k,L). Suppose also that LAY = 10 so that the 
dimensions of the matrix are 10 x 10. The initial values 
of I(T) su that oie pa the matrix which will be used 


i ty ip 


to peterate the permutations are shown in Table 8. 
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Table 8. Initialization of Matrix N(k,L) 
and Vector J(I). 


Column > == 7718. 9) 10 J(7) = 7 Original Partition 
Row ¥ | 
J(8) = 8 2133: 1) 
7 2€G)|0¥3 3 
8 eres 3 J(9) = 9 
9 2 ates 3 
10 AN ee 3 


The value of J(9) is increased by one to 10 so that 
the integers are to be interchanged in positions N(9,9) and 
N(9,0(9)), that is, in positions N(9,9) and N(9,10). They 
are equal therefore no exchange is made. The value J(8) is 
now increased by one to 9 and J(9) is returned to its origi- 
nal value of 9. Integers in positions N(8,8) and N(8,J(8)) 
are to be interchanged. This operation is performed and the 
resulting permutation is written into rows 9 and 10 as shown 


in Table 9. Row 8 remains unaltered. J(9) is now increased 


Table 9, Generation of Sequence (2) 
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by one to 10 and the integers in positions N(9,9) and 
N(9,J(9)) are to be interchanged. The operation is per- 
formed and the resulting permutation is written into row 10 
as shown in Table 10. Row 9 remains unchanged. J(8) is 

now increased by 1 to 10 and J(9) is returned to its' origi- 
nal value. Integers in positions N(8,8) and N(8,J(8)) are 


to be interchanged. No interchange is made due to the second 


Table 10. Generation of Sequence (3) 


——. 


Column > PRBO GR eto dG7') Saee Wia Permutation Generated 
Row ¥ a) 
J(8) = 9 23°30 14639) 
7 Zyl 33 3 
8 Catal Lane 3 J(9) = 10 
9 Voges Foe | 3 
10 2i7 3) 3 ] 


criterion given previously for the generation of distinct 
permutations. If the interchange were made, the result 
would be identical to sequence (3). Since J(8) cannot be 
increased beyond 10, it is returned to its' initial value 

of 8 and J(7) is increased by one to 8. Integers in posi- 
tions N(7,7) and N(7,J(7)) are to be interchanged. Once 
this is done all distinct permutations involving integers 

in colums 8, 9 and 10 are generated as previously shown. 

The value of J(7) is then increased by one and the procedure 
is repeated. This continues until all permutations involv- 


ing J(7) = 10 have been generated. Following this example, 
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the operation of the algorithm should be clear. 


Ae See COSC Comparison Between Hron's Algorithm and the 
ermute artition Approac 


The total amount of vector space necessary. for the 
generation of the kinematic code using the permuted parti- 
tion approach is given by (LAY)* + 3(LAY). These values 
are tabulated and compared with the space needed for Hron's 
matrix N(I,L) in Table 11. As the number of layers is 


increased it can be seen that the saving in space becomes 


BaD ewe Comparison of Vector Space Requirements 
for the Generation of the Kinematic Code. 


LAY N(I,L) (LAY) + 3(LAY) 
2 4 10 
4 24 28 
6 120 54 
8 560 88 
10 2,520 130 
12 11,088 180 HMS = LAY 
14 48.048 238 
16 205,920 304 
18 875,160 378 
20 3,695,120 460 
22 Nes 1o. 504 550 
24 64,899,744 648 
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enormous using the permuted partition approach. The reduced 
amount of space should result in a more economical means of 
generating synthetic seismograms. In order to test this, 


two computer programs were written to generate synthetic 
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seismograms at vertical incidence. They were identical in 
all respects except in their method of generating the kine- 
matic code. One program used the algorithm of Hron (1972), 
the other made use of the permuted partition approach. For 
this test the maximum allowable number of half-setments (HMS) 
was set equal to LAY. The resulting cost breakdown is shown 
in Table 12. 

For models of twelve layers or less the space necess- 
ary for the generation of the kinematic code in both programs 
is considerably less than the necessary peripheral space; 
that is the space required for the model parameters, spike 
synthetic, wavelet storage, and the plotting routines. 
Therefore in the case of the nine and twelve layer models 
the difference in storage costs between the two programs is 
rather small. However, in the case of the fourteen and six- 
teen layer models, the required space for the matrix N(I,L) 
has grown very rapidly beyond the peripheral storage require- 
ments. In addition, the difference between the values 
(LAY)® + 3(LAY) and (I)(L) for LAY = 16 is much larger than 
the corresponding difference for LAY = 9. The result is a 
reduction in total computing costs by approximately one- 
half when using the permuted partition approach for sixteen 


layers and HMS = LAY. 


2.6 Restricted Partitioning 


Despite the reduction in cost obtained by using the 


method of permuted partitions the generation of the synthetic 
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seismogram for sixteen layers with HMS = LAY was still 
relatively expensive. It would be useful therefore to 
examine the effect of further restricting the number of 
phases that contribute to the final result. 

It should be noted that the ray selection criteria 
of Hron (1972) used in the cost comparison studies are 
themselves a restriction, as is setting the value of HMS 
equal to the number of layers (LAY). The result is an 
incomplete or partial ray expansion, where only those rays 
that are thought to contribute large amplitudes are included. 
M. Hron (1973) studied the problem of ray series convergence 
one detail and was unable to find any theoretical 
criteria by which the ray series could be terminated with 
assurance that the remainder would be negligible. The 
result of this is that any partial ray expansion must be 
the result of intuition as to which group of rays contribute 
the maximum amplitude. With this in mind further restric- 
tions on the partial ray expansion of Hron (1972) are 
examined. 

In order to reduce the number of phases, restrictions 
on the number of parts allowed in a given partition were 
imposed in four different levels of severity. These are 
shown in Figure 10. The partitions to the right of the 
stepped line in each of Figure 10(a), (b), (c) and (d) are 
used in the computation while the partitions to the left 
are discarded. This approach tends to eliminate phases 


with a large number of reflections by disallowing the distri- 
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(c) Severity Level 3 
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(b) Severity Level 2 
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(d) Severity Level 4 
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bution of a large number of half-segments into a small 
number of layers. The restriction on the number of reflec- 
tions could of course be applied directly when generating 
the dynamic portion of the ray code, however it would then 
be necessary to produce the complete set of partitions. 

In the method described the suspect partitions are elimin- 
ated directly by simply not generating them resulting in a 
greater saving of CPU time. 

The restricted partitioning was first applied to a 
sixteen layer model, (Figure 11), which approximates the 
geologic structure near Edmonton. The true structure of 
the area was obtained from a digitized sonic log of well 
CS Big Hay Lake located at 13-30-48-22W4. The log provided 
data for the interval between 670 and 6,030 feet. The 
problem of how to best approximate a complex velocity struc- 
ture such as that obtained from the digitized sonic log by 
a model with a small number of layers is not dealt with here. 

In Figure 12 the spike synthetics obtained using the 
sixteen layer model are shown for the unrestricted and four 
different levels of restricted partitions. The total number 
of different phases represented in each of these spike 
sequences along with the costs involved in producing the 
final synthetic seismograms are presented in Table 13. 

There is a noticeable change in the character of the spike 
sequences in progressing from the unrestricted case to 
Severity Level 4. In particular there is a large reduction 


in the number of phases arriving at times greater than 1.105 
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Table 13. Cost Comparison of Four Levels of 
Restriction in the Generation of 
the 16-Layer Synthetic Seismogram. 


NS 
a 


Total Number CPU Time CPU Samant 


of Phases (sec. ) (pg-min. 
Unrestricted 48,760,342 254.8 292.0 
Severity Level 1 87 ,668 17.28 19.4 
Severity Level 2 16,343 9.6 10.0 
Severity Level 3 Da) 720 ne 
Severity Level 4 24] D9 ey) 


seconds which marks the end of the velocity structure mea- 
sured in two-way travel time. Hence the restricted parti- 
tions have been successful in terms of reducing the number 
of late arriving phases with a corresponding drastic reduc- 
tion in overall computing costs. However, there has also 

been a reduction in the number of phases arriving at times 
of less than 1.105 seconds and the question remains as to 

how severe this reduction has been. In order to obtain a 

quantitative answer to this question a relative percentage 
difference was calculated for each .2 second time interval. 


That is, for each interval: 


Average Relative dd y i Te 
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where U. is the amplitude value of the spike in the un- 
restricted sequence at time 1/1000 and Rs is the amplitude 
value of the spike in a restricted sequence at the identical 
time; Umax is the amplitude value of the largest spike from 
the unrestricted sequence in the time interval, and n is the 
number of non-zero values of U. within the interval. If the 
value U. is found to be equal to zero the operation 

|U. - Rel/[Uyay| is not performed. The reason for using 

the expression |U, - R.|/|Uyaylin (16) rather than the usual 
percentage difference given by JU. - R.|/|U,| is that the 
latter does not yield a good measure of the dissimilarity 
between the two final synthetic seismograms. This is because 
the absence of very small spikes in a region characterized 
by spikes of much larger amplitude, has only minor effects 
on the final trace, yet the values |U. - R.|/|U.| could be 
very large giving a false indication of poor agreement 
between the two seismograms. The measure in (16) was app- 
lied to all four severity levels in Figure 12 and the 
results are shown in Table 14. The bracketed quantities 
give the maximum relative percentage differences within 

each time interval. Some of these have high values which 
indicate that there are spikes of significant amplitude 
which have been excluded by the restriction process. The 
values in the last column of Table 14 give a quantitative 
measure of the spike sequence degeneration as seen in 

Figure 12. To investigate the dependence of the averages 


over 1.105 seconds on the length of the time interval chosen, 
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Table 14 was recalculated using a time interval of .1 
seconds. The changes noted were not significant. 

The final synthetic seismograms were produced by 
convolving a wavelet with the spike sequences in Figure 12. 
These are shown along with the sonigram in Figure 13. Al- 
though the general features are present throughout, the 
Character changes in going from the unrestricted wave train 


to Severity Level 4. 


2.7 Synthetic Seismograms with Geometric Spreading 


The spike sequences in Figure 12 were generated 
without taking into account the factor introduced into 
the amplitude as a result of geometric spreading of the 
wavefront. One advantage of the ray tracing method is that 
inclusion of this factor is a relatively simple matter. 
Synthetic seismograms based on a 14-layer model were gene- 
rated both with and without geometric spreading in order to 
observe any changes in character. The spike sequences were 
generated letting HMS = LAY and the wavelet used in generat- 
ing the final traces was identical to the one shown in 
Figure 12. The spike synthetics and convolved traces that 
resulted from unrestricted partitioning are shown in Figure 
14 along with the sonigram. 

It is clear that if the spreading correction is app- 
lied to the spike sequence in Figure 14(a), the changes in 
amplitude of the phases will not occur in a uniform or 


predictable manner when viewed in time. This is confirmed 
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by noting the large change in relative amplitudes of the 
group of arrivals at .2 and .5 seconds in Figure 14(b) as 
compared with the same arrivals in Figure 14(a). In the 
remainder of the spike sequence however, the differences: 
are rather subtle. 

The measure of dissimilarity used in the last section 
was applied to the pestis of restricted partitioning with 
geometric spreading using the 14-layer model. The process 
was then repeated without the spreading correction. Table 


15 displays the results. 


Table 15. The Average Relative Percentage 
Difference With and Without 
Geometric Spreading. 


Average Relative Percentage Differences 
Over 1.105 Seconds 


Without Spreading With Spreading 
Severity Level 1 wi A9 306 
Severity Level 2 4.7 Aa 
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CHAPTER III 


THE RAY METHOD OF Nth ORDER MULTIPLES 


3.1 Basic Concepts 

The last two chapters were centered around the ray 
method of Hron (1972). In this chapter the classification 
scheme commonly used in industry is described. This 
classification is based on the number of reflections and 


rh order 


the groups of ray paths which result are called n 
multiples. The order of the multiple is defined as the 
number of downward reflections suffered by the ray in 
travelling from the source to receiver; that is, the number 
of times an upgoing ray segment is transformed into a down 


going ray segment by a reflection from an interface. Figure 


15 may help to illustrate the concept. The ray paths shown 


Figure 15. First and Second Order Multiples. 
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in Figure 15(a) and (b) are both first order multiples 
while those in (c) and (d) are second order multiples. 
Although the method described in Chapter I and the 
walt Order multiple approach are based on entirely different 
criteria, the two sets containing the phases for each scheme 
using a j-layer velocity model, will in general have a 
fairly large intersection. The degree to which this hold 
true depends on the specification of the value HMS and the 
highest order of multiples to be used. It is entirely 
possible for one set to exist as a subset of the other. 
For example by specifying HMS = 2 * LAY, the first order 
multiples are automatically included using the scheme 


described in the first chapter. Letting HMS = 4 * LAY 


includes first and second order multiples. 


3.2 The Coding Scheme and Method of Generation 


The code for a particular ray path using the approach 
of nth order multiples consists of 2n + 1 integers where 
n represents the order of the multiple. Each integer m(i) 
of the code represents the interface from which the yee 
reflection occurs. Some typical ray paths and their codes 
have already been seen in Figure 15. For first order 
multiples (Figure 15(a) and (b)) the range of values for 
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The range of values for the rg 


order multiples may be 
specified in a similar manner. 

In order to generate the codes, the sequence of 
inequalities in (17) were written into nested "do-loops". 
Using this method, vector space is needed only for storing 
the. code currently being generated. As in the scheme 
described in Chapters I and II, each code completely speci- 
fies the kinematic and dynamic properties of a single phase 
travelling from the source to receiver. However in the nth 
order multiple approach at vertical incidence using the 
code generation scheme described, it is not necessary to 
determine the amplitude of each of these phases separately. 
In the computer program written to generate synthetic seismo- 
grams using this approach, the amplitude of the phase 
currently being generated is obtained from the amplitude 
of the previous phase. To understand how this is possible 
consider the two phases drawn in Figure 16. According to 
the scheme used, the phase shown in Figure 16(b) is gene- 
rated immediately following the phase shown in (a). To 
obtain the amplitude of the phase in (b) it is necessary 


to divide the amplitude of the phase in (a) by the reflec- 


tion coefficient at interface 2, multiply by the downward 
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Figure 16. Two Sequential First Order Multiples. 


(a) (b) 


and upward transmission coefficients at interface 2 and 
finanay multiply by the reflection coefficient at interface 
3. Since calculating the amplitude in (b) would require 
eight operations if performed independently of (a), a sav- 
ing of four operations has been realized. The savings 
increase as the higher order multiples are included. 
Unfortunately this saving may be overshadowed by the fact 
that the concept of dynamic equivalence has not been applied 


to this scheme. 


3.3 Comparison of the Method of Nt" Order Multiples 
and the Ray Method of Hron (1972) 

A comparison was made between the two different ray 
classification schemes using the 16-layer model shown in 
Figure 11. The spike sequences for primaries plus up to 
second order multiples are shown in Figure 17 together with 


the Severity Level 3 spike synthetic. The average percent- 
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age differences indicated are relative to the Unrestricted 
sequence in Figure 12. These values for each .2 second 
interval up to the end of the velocity structure at 1.105 
seconds, are shown in Table 16. Bracketed quantities in 
Table 16 give the maximum relative a aiariaians differences 
encountered in the time intervals. 

It is interesting to note that the number of phases 
in the primaries plus first and second order multiples spike 
sequence is almost two orders of magnitude larger than the 
number of phases in the Severity Level 3 sequence. Despite 
this, the Severity Level 3 sequence has a smaller relative 
percentage difference over 1.105 eeeonde The reason is 
that most of the contributions of the primaries plus first 
and second order multiples occur at times greater than Pate 
seconds. 

The wavelet shown in Figure 13 was convolved with the 
spike sequences in Figure 17 to yield the synthetic seismo- 
grams displayed in Figure 18. The costs of generating 
these synthetics is shown in Table 17. The effective cost 
per phase in generating the Severity Level 3 sequence is 


° seconds per 


SOx 107? seconds compared with 6.1 x 10° 
phase for the primaries plus first and second order multiple 
sequence. The large difference occurred despite the fact 
that the generation of the primaries plus first and second 
order multiples sequence was done without making use of 


dynamic equivalence. It might be speculated that if restric- 


tions were imposed on the primaries plus first and second 
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Figure 18. Synthetic Seismograms for Severity Level 3 and 


Primaries Plus up to Second order Multiples 
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Table 17. Cost Comparison of Generating 
Primaries Plus Up to Second 
Order Multiples for the 16- 
Layer Model. 


Se 
SS 


Total Number CPU Time CPU ieeteit 


Of Phases (sec. ) (pg-min. 
Severity Level 3 25335 . J0 1 ey 
Primaries Plus Ist 
and 2nd Order 164,576 10% 0 10.9 
Primaries Plus Ist 
Order Te ti2 be 6.0 
Primaries Only 16 1.9 Dc 


order multiples generating scheme in such a way that more 
of the phases were chosen to lie within the 1.105 second 


th order 


time window, then a restricted version of the n 
multiple approach might be more economical than the scheme 
described in Chapters I and II. Unfortunately no investiga- 
tions were made in this area. It may be of interest to note 
that the effective cost per phase in generating the 


Unrestricted sequence in Figure 12 using the concept of 


dynamic equivalence was 5.2 x 107° seconds. 


3.4 Surface Order Multiples 


yth order "surface" multiples occur as a subset of 
nee order multiples and are defined as those multiples in 
which m(i) = 0 for i even. For examples refer to Figure 


15(b) and (d). Surface multiples have often been used as 


0s le t steet 


; ae 
ey wrod? usd "eee 59, wee Areker p 
mtm=ug) (384) ie: 2d tO” 


1 Net | abe t Teves wl rrevee 


Ve | ted aut a 
@.01 | oor “avd bar my Pa 
| “ger auld estapehet " 
: v8by0 


0.8 Weta Ae ) 
| EU col Bt mei a ein soeaainns | Oe 
sSrEs ie ae ett ai ™ aay 


stom tent New 5: “Adu iif pinarae) piri tovenog 2otqts um sabre 
bnoose Zor. ‘ wits ‘ahhh avi oF ne2ono evow. aseedg add to. 
yabr0 3 n ag to nobevay. bagatit 20% 5 ners ewobn hw emts— ve 
suena end nse aatmengse stom od Siptm damovaas ofah3stom t, ae 
qbetizevat: on vied onuayotny TT bing, a | vat qaild at bedinsesb ie 
agon ot taenadnt to od ‘6m + som atad of abou atom. enots ae 
ong ectiyevensp: ve g2eng en, $200 svttoeyte ins todd ove 

to Sqeane> aad ‘pniew St equip nt someupad betotrtesynh. i 

Lae ee a . an0982) mS “Or. Xue To Lesa , 


16 bn zataraton aire Oy 


“= a er. cee 


3.6) 


an approximation to the complete set in designing multiple 
Suppression algorithms, e.g., Watson (1965). To observe 
how well the approximation holds a comparison was performed 


between nth order multiples and nth 


order surface multiples 
up to the second order using the 16-layer model. Figure 19 
shows the spike sequence and Figure 20 the final synthetic 
seismograms. The average percentage differences indicated 
in Figure 19 are relative to the Unrestricted spike sequence 
in Figure 12. The wavelet used to generate the synthetics 
was identical to the one shown in Figure 13. Agreement 
between the surface multiples and the complete set up to 


the second order is seen to be quite good in the case of the 


model used. 


3.5 The Complete Wave Solution for the 16-Layer Model 


To this point in discussing the various restrictions 
and schemes used in generating spike sequences, comparisons 
were generally made with the Unrestricted sequence in Figure 
12. The reason for this is that this sequence represented 
the effect of the largest number of phases and was assumed 
to be the best approximation to the complete solution 
obtained by distributing up to an infinite number of half 
segments into a maximum of LMX layers. However if the 
Unrestricted sequence does not adequately represent the 
complete solution, then the comparisons made are of little 
Significance.p bor this reason a spike sequence representing 


the complete wave solution was obtained and compared with 
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Figure 20. Comparison of ayn Ge Seismograms for ee Order 
h order Surface Multiples 
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the Unrestricted sequence in Figure 12. These are shown 
in Figure 21 together with the synthetic seismograms. As 
can be seen, within the time window of 1.105 seconds the 
agreement is very good. This relieves any doubt which may 
have existed as to how meaningful previous comparisons have 


been. 


3.6 Conclusions 

The new approach developed for the generation of the 
ray code using the classification scheme of Hron (1972) 
proved to be successful in reducing computer storage costs. 
For the sixteen layer model with HMS = LAY, the reduction 
in total cost for generating the synthetic seismogram was 
Slightly greater than one-half. Restricted partitioning 
further reduced costs by generally suppressing those contri- 
butions outside the time range of the velocity structure. 
The Severity Level 1 restriction resulted in a synthetic 
seismogram at approximately 1/20th the cost required to pro- 
duce the Unrestricted sequence using the code generation 
scheme described in Chapter I. In addition the Severity 
Level 1 sequence retained most of the character observed 
in the Unrestricted sequence within the time window of 1.105 
seconds which marks the end of the velocity structure in two- 
way travel time. 

Inclusion of the geometric spreading correction in 
the amplitude computation was found to significantly affect 


the character of the synthetic seismograms, particularly 


A 


van atau $du0d Yrs apie att : 
avon anoetysqmos eug ter: iwtonraneh ait 


A = | 
aft to notisevensp ait oF beyolsvab iononaas wert edt: by WE 
(syer) nov +0 omantde | gels i tibeing ort pniew ahaa yoy 


,2t2oo Bus yOtya ver wadtoD patyuber nt lutadeooue od ae bevora . 7 
fotsaybay ert qY Rus an Tas rut Seber wave (mesante oft N04 iva 
25W davpometsd af bddgade: etd pnitaronsy. 16% F209 fare nt un 
ontnot¢iis+eg hatotyt2er Thad NO, neds Tavsarg ylanette ap 
-~trtnod s2e0dd diidinlaliaiaalaen ud et209 | abies yeti tyu? av es 
srvtouyse ytisolsy. ent 10 epasy gms oid abFedue anottud | ie 
ottadtaye 6 nf botfuest nottotrtes I fevet yttreve2 onT ; ois 
-07q od bevtupet teo9 ana asOS\T viatomtxo1ggs Je merpometez | * 
noftevenee abs’ ott ghtelr nuKueen hetoriteer nv. ott eoab; ) a 
ysineve2 odd vor dtbbs al y Det qonto wt badtroesd amatioa) a | : 
bavrsedo vetsersits: oat to, Jeom bantstey saneupse | faved : 
201.1 to wobatw ants 849 Abidtw isonsupes: Botatat sora aft mt u . 
Bid nt oxut ound, xt 199 ti aa lee ond, advem dotnw detache ! 
en emit Povend Yow Pe | 
p ap oe a2 to aoreutont uy he r 


sige moktvingnos ‘bud Figs ek 


59 


J3S-3WIL 

Ov°2 02°22 0o°2 og°t o9°1t Or’ tl 02° tT oo°T 0e°o 09°0 ov°o 02°0 00°0, 
= 

°np 

NOILNIOS 31L37dNOD z 

of, 

om 

(=) om 

oOo 

Mm 
G3i21¥LS3YNN 2 


NOILMIOS 3L371dROD 


A CATAVAVAVAVATAUAUITANAVAVAVACAUAVAUA A cacaacacaemcacca haere’ 


Q31L0INLS3SYNN 


Tepow zoAeT-9T By AOF SuOTANTOS 
SAPM SReTdUOD pue paezoTAZSeAUN sy FO uoSTAedUOD “Tz sanbty 


as ee nen 


“4 + ae us ‘i i on Ws ble ine i 
. f] . i ” 
- 4 7) : » 


60 


the relative amplitude of different groups of arrivals. 


In comparing the method of nth 


order multiples with 
the classification scheme of Hron (1972), it was found that 
for the sixteen layer model, the synthetics including pri- 
maries plus first and second order multiples compared 
closely with the Severity Level 3 sequence within the time 


window of 1.105 seconds. However the nth 


order multiple 
approach seems to be less selective in the rays that it 
generates in the sense that a large majority of these rays 
fall outside the time window. 


nth order surface multiples were found to be a good 


approximation to the complete set of nih 


order multiples 
for the model used in the comparison. 

The Unrestricted sequence for the sixteen layer model 
was seen to correspond very well with the complete wave 
solution thereby justifying previous comparisons. 

Some time after completing the work on the permuted 
partition approach it was discovered that a similar techni- 
que had been successfully applied by Vered and Ben-Menahem 


(1974). However, they do not appear to make use of 


restricted partitioning to reduce the number of undesir- 


able phases. 
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APPENDIX A 


NUMBER OF ARRANGEMENTS OF n, IDENTICAL BALLS 


2 
INTO ny POCKETS, SOME OF WHICH CAN BE EMPTY 


(FROM HRON (1972)) 


Consider an arbitrary distribution of No identical 
balls into n, 2 1 pockets separated by (n, - 1) inner walls. 
Each such distribution can be visualized fa a sequence of 
(n, - 1 + n,) items represented by n, balls and (n, - 1) 


inner boundaries, as illustrated in Figure A.1 for n, = Ar 


Noa 5. From this diagram, we can see that three inner 
Figure A.1. One Distribution of 5 Balls 
into 4 Pockets. 
I IL Il 
Pe ea 
2 3 4 5 6 7 


walls occupy, in the sequence of eight elements, the fourth, 
fifth, and seventh positions so that there are 3 balls in 
the first pocket, the second pocket is left empty, and there 


is 1 ball in pockets III and IV. 


Because individual distributions differ from each 
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other only in selection of (n, - 1) positions occupied by 
inner boundaries, the number, My» of distributions of No 
balls into ny pockets, some of which can be left empty, 
is equal to the number of all possible selections of 

(n, - 1) positions from the total number (es 1b no). 


Obviously, this is equal to 


n,+n,-1 n,+n,-1 (ny, © Nase l jl 
Dae ere 1 z 
M,(n,,n,) = C c = A. 
ie pa for 
There C, _y is the number of all combinations of the 
1 


(n, - 1)th class from (n, sind 1) elements. 
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APPENDIX B 


NUMBER OF ARRANGEMENTS OF n, IDENTICAL BALLS 


2 
INTO ny POCKETS, NONE OF WHICH CAN BE EMPTY 


(FROM HRON (1972)) 


Distributions of this kind can be performed only if 
ny 2 M,. Then the problem can be transformed easily to the 
problem resolved in the previous section if we separate 


from n, balls ny so-called "fixed balls" by distributing 


2 
them into ny pockets, none of them leaving empty. Then the 
number M, has to be equal to the number of arrangements of 


the remaining (no - ny) "movable balls" into ny pockets 


provided that some of them may not contain any movable ball. 


Using Equation A.1.1 we get immediately 


no-l (n5--. 1)! 


ba - 2h eee £ 
M(n,.n) a My (ny sn, Fi ny) r Cn -1 = (Hy = )ith, = naif = i no = il, I Bei 


where M,(n,.n5) is the number of all arrangements of No 


identical balls into ny pockets leaving none of them empty. 
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APPENDIX ¢ 
COMPUTER PROGRAMS 


C.1 N-th Order Multiples Approach 

3 FORMAT (//,T10,*DOWNWARD REFLECTION COEFFICIENT!,//) 
12 FORMAT (7F10.5) 

4 FORMAT (//, 710, "DOWNWARD TRANSMISSION COEFFICIENT',//) 
13 FORMAT (7F10.5) 

5 FORMAT (//,T10,"UPWARD REFLECTION COEFFICIENT? 7/7) 
14 FORMAT (7F10.5) 

6 FORMAT(//,T10,"UPWARD TRANSMISSION COEFFICIENT!,//) 
15 FORMAT (7F10.5) 

7 FORMAT (*TIME=',F10.5) 

8 FORMAT ("AMP=",E£15.4) 


11 FORMAT (//,T28,*PRIMARY REFLECTIONS',//) 
16 FORMAT (//,"GROUP SYMBOL',T30,15) 

17. FORMAT (T40,15) 

an FORMAT (4F10.5) 

22 FORMAT (5F10.2) 

23 FORMAT(F10.2) 


504 FORMAT(/,'TOTAL NUMBER OF PHASES',1T30,110) 


DIMENSION THICK (16) 
DIMENSION VEL(17) 
DIMENSION DEPTH(16) 
DIMENSION TIME (16) 
VALUE CONTAINS THE WAVELET PLUS TRAILING 
ZEROS. THE DIMENSION OF VALUE MUST BE GREATER THAN OR 
EQUAL TO THE LENGTH OF THE CONVOLUTION OF THE WAVELET 
WITH THE SPIKE SYNTHETIC, WITH THE ADDED CONDITION THAT 
THE DIMENSION MUST EQUAL 2**N WHERE N IS SOME POSITIVE 
INTERGER. 2**5=32 2**6=64 2**7=128 2**8=256 2**9=512 
2**10=1024 2**11=2048 2**12=4090 2**13=8192 2**14=16384 
COMPLEX VALUE (8192) 
P50 CONTAINS THE SPIKE SYNTHETIC (OUTPUT FROM RDOFF) 
DIMENSION P50(SAME DIMENSION AS VALUE) 
COMPLEX P50(8192) | 
DIMENSION P60(4602) 
DIMENSION BNUIS(2 * DIMENSION OF VALUE) 
DIMENSION BNUIS (16384) 
FQUIV ALENCE (VALUE, BNUIS) 
C CONTAINS THE Y OR AMPLITUDE VALUES USED IN PLOTING THE 
SYNTHETIC SEISMOGRAM 
DIMENSION C(LENGTH OF THE CONVOLUTION OF THE SPIKE SYNTHETIC 
WITH THE WAVELET + 2) THE SCALING PARAMETERS ARE STORED IN 
THE LAST TWO POSITIONS 
DIMENSION C (4602) 
C201 CONTAINS THE X OR TIME VALUES USED IN PLOTTING THE SS 
DIMENSION C201(SAME AS DIMENSION OF C) 
DIMENSION C201(4602) 
IN NLOGN THE DATA LENGTH MUST BE EQUAL TO 2**N WHERE N IS 
SOME POSITIVE INTERGER. DIMENSION MMM(LARGEST VALUE OF N TO 


BE PROCESSED) 
DIMENSION MMM(13) 


GTIME(I) REPRESENTS THE VERTICAL TRAVEL TIME FROM THE SURFACE 
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TO THE BOTTOM OF THE ITH LAYER 
DIMENSION GTIME(131) 
COMPLEX B,CMPLX 


DO 45 I=1,8192 
P50 (L) =CMPLX(0.0,0.0) 
45 CONTINUE 


LAY=16 

JJSO=LAY+1 
INPUT DATA 

READ(7,23) (DEPTH (I) ,T=1,LAY) 

READ(8,23) (VEL(I) ,I=1,JJ50) 
IF THE GEOMETRICAL SPREADING FACTOR IN THE AMPLITUDE IS NOT 
DESIRED SPECIFY TIG=0. ANY OTHER VALUE FOR TIIG 
WILL INVOKE THE SPREADING CALCULATION. BE CERTAIN 
TO CHANGE THE VERTICAL SCALING PARAMETERS IN THE 
PLOTTING ROUTINE WHEN RUNNING THIS PROGRAM 
WITHOUT GEOMETRICAL SPREADING 

IIG=0 . 
IF FIRST ORDER MULTIPLES ARE NOT DESIRED 
SPECIFY. .LLF=0 

ITF=1 
IF SECOND ORDER MULTIPLES ARE NOT DESTRED 
SPECIFY IIS=0 

IIs=1 
TO SUPPRESS PRINTOUT OF REFLECTION AND TRANSMISSION. 
COEFFICIENTS SPECIFY ISUPCO=0. ALL OTHER VALUES 
WILL RESULT IN A PRINTOUT 

ISUPCO=0 
TO SUPPRESS PRINTOUT OF AMPLITUDES AND TIMES 
FOR DIRECT REFLECTIONS SPECIFY ISUPP=0. ALL 
OTHER VALUES WILL RESULT IN A PRINTOUT 

ISUPP=0 
TO SUPPRESS PRINTOUT OF AMPLITUDE AND TIMES FOR 
FIRST ORDER MULTIPLES SPECIFY ISUPFO=0. ALL OTHER 
VALUFS WILL RESULT IN A PRINTOUT 

ISUPFO=0 
TO SUPPRESS PRINTOUT OF AMPLITUD AND TIMES FOR 
SECOND ORDER MULTIPLES SPECIFY ISUPSO=0. ALL OTHER 
VALUES WILL RESULT IN A PRINTOUT 

ISUPSO=0 
TO SUPPRESS PLOT SPECIFY ISUPLOT=0. ALL OTHER 
VALUES WILL RESULT IN A PLOT OF THE SYNTHETIC 

TSUPLT=1 

THICK (1) =DEPTH (1) 

M=LAY-1 

M1=LAY-2 


DOea Venus 2 puAy 
THICK (I) = DEPTH (I-DEPTH (T-1) 
10 CONTINUE 


po 20 I=1,LAY 
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TIME(I) = THICK(I) /VEL(I) 
CONTINUE 


CALCULATE THE REFLECTION AND TRANSMISSION COEFFICIENTS 
DIMENSION DREFCO (16) 
DIMENSION DTRACO(16) 
DIMENSION UREFCO(16) 
DIMENSION UTRACO(16) 


DO 25 I=1,LAY 
DREFCO (1)=(VEL(I-VEL(I+1))/(VEL(I) # VEL(I+1)) 
CONTINUE 


IF(ISUPCO .FQ. 0) GO TO 4y 
WRITE (6,3) 
WRITE(6,21) (DREFCO(I),1=1,LAY) 


DO 30 I=1,LAY 
DTRACO(I)=1+DREFCO (I) 
CONTINUE 


IF (ISUPCO.FQ.0) GO TO 46 

WRITE (6,4) 

WRITE(6,21) (DTRACO(I) ,I=1,LAY) 
UREFCO(1)=1.0 


DO 35 I=2,LAY 

UREFCO (I) = (VEL (I-VEL (I-1)) / (VEL (1-1) #VEL(I)) 
UREFCO INDICATES REFL. FR. TOP OF ITH LAYER 

CONTINUF 


IF (ISUPCO.EQ.0) GO TO 47 
WRITE (6,5) 
WRITE(6,21) (UREFCO(I) ,I=1,LAY) 


DO 40 I=2,LAY 
UTRACO(I) =1+UREFCO (I) 
CONTINUE 


UTRACO(1)=1.0 

IF (ISUPCO.EQ.0) GO TO 48 

WRITE (6,6) 

WRITE (6,21) (UTRACO(I) ,I=1,LAY) 
IAFG=0 

NOW CALCULATE TIME AND AMPLITUDE FOR DIRECT REFLECTIONS 
AMP=1.0 

IF (ISUPP.EQ.0) GO TO 43 

WRITE (6,11) 

IND=1 

TIM=2¥*TIME (1) 

IP(IIG .EQ. 0) GO TO 41 

AM P= AMP* DREFCO (1) / (2* DEPTH (1) ) 


GO TO 42 
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44 AM P= AMP*DREFCO(1) 

42  IF(ISUPP.EQ.0) Go TO 44 
WRITE (6,16) IND 
WRITE(6,17) IND 
WRITE(6,8) AMP 
WRITE(6,7) TIM 

44 B=CMPLX (AMP, TIM) 

CALL RDOFF(B,P50,IAFG) 
GTIME (1) =TIME(1) 


DO 80 I=2,LAY 
AMP=1.0 
T1=I-1 


pow6S J=1, 11 
AMP=AMP*DTRACO(J) 
65 CONTINUE 


DOV, T 
AMP=AMP*UTRACO (J) 
75 CONTINUE 


TIM=TIM+ 2*TIME (I) 

IND=IND+1 

IF(IIG .£Q. 0) GO TO 76 

PL=0. 

po’'78} d=1,1 

PL=PL+2* THICK (I) /VEL(1) 
78 CONTINUE 

AMP=AMP*DREFCO (I) /PL 


GO TO 77 


76 AMP=AMP*DREFCO(T) 
77. = GIME( 1) =TIM/2 
B=CMPLX (AMP, TIM) 
IF (ISUPP.EQ.0) GO TO 86 
WRITE(6,16) IND 
WRITE(6,17) IND 
WRITE (6,8) AMP 
WRITE(6,7) TIM 
86 CALL RDOFF(B,?50,IAFG) 
80 CONTINUE 


IF(IIF.EQ.0) GO TO 81 
CALL FOM(IND,LAY,TIME, GTIME,DREFCO,THICK, DEPTH,IIG, 


CDTRACO, UREFCO, UTRACO, ISUPFO,B, P50, TAFG) ~ 


81 IP(IIs.£0.0) GO TO 82 
CALL SOM (IND, LAY, TIME,GTIME, DREFCO, THICK, DEPTH, IIG, 


CDTRACO, UREFCO, UTRACO, ISUPSO,B, P50, LAFG) 
GENERATE THE TIME AXIS FOR PLOTTING 
IF (ISUPLT.EQ.0) GO TO 2004 
Ape c20TC1)=0.0 
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C INSERT DO LOOP PARAMATER. I=2,LENGTH OF TRACE 
C IN SECONDS X 1000 


DO 2001 T=1,4600 
C201(1I) =I/1000. 
2001 CONTINUE 


C SPECIFY THE KLAUDER WAVELET PARAMETERS 
Cc CALL KLAUD(7.,12.,4., VALUE) 
CALL WAVLT1(VALUE) 
DO 83 I=1,4602 
P60 (I) =REAL (P50 (I) } 
83 CONTINUE 
WRITE(11) P60 
C IN THE NEXT TWO STATEMENTS THE FIRST PARAMETER OF 
C NLOGN MUST EQUAL THE DIMENSION OF MMMQ ) 
CALL NLOGN(13,VALUE,-1.0) 
CALL NLOGN(13,250,-1.0) 


DO 2002 I=1,8192 
VALUE (TI) =VALUE(T) *P50(T) 
2002 CONTINUE 


C INSERT 1ST NLOGN PARAMETER 

CALL NLOGN(13,VALUE,1.0) 
C INSERT DO LOOP PARAMETER, I=1,LENGTH OF TRACE 
C IN SECONDS X 1000 


DO 2003 T=1,4600 
C (I) =BNUIS (2*I-1) 
2003 CONTINUE 


WRITE(6,504) IAFG 

CALL PSY (C,C201, P60) 
2004 STOP 

END 


C *FOM* CALCULATES THE TIME ARRIVAL TIM AND AMPLITUDE 
C OF THE FIRST ORDER MULTIPLES 
SUBROUTINE FOM(IND, LAY, TIME, GTIME, DREFCO, THICK, DEPTH, IIG, 
CDTRACO, UREFCO, UTRACO, ISUPFO,B, P50,IAFG) 
DIMENSION GTIME(131) 
COMPLEX B,CMPLX 
COMPLEX P50(8192) 
DIMENSION DREFCO (131) 
DIMENSION DTRACO (131) 
DIMENSION UREFCO (131) 
DIMENSION UTRACO (131) 
DIMENSION TIME(131) 
DIMENSION THICK (131) 
DIMENSION DEPTH (131) 
100 FORMAT (//,T28,"FIRST ORDER MULTIPLES',//) 
101. FORMAT(//,'GROUP NUMBER',1T30,15) 
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FORMAT (*+",T40, 314) 
FORMAT ("AMP =", §£15.4) 


| FORMAT ("TIME =", F10.5) 


FORMAT (F 10.5) 

AMP3=1.0 

IF (ISUPFO.EQ.0) GO TO 110 
WRITE (6, 100) 


DO 200 IA=1,LAY 
PTIME=GTIME(IA) 

GSP=0. 

DO 210 I=1,IA 
GSP=GSP+THICK (I) *VEL(I) /VEL (1) 


' CONTINUE 


AMPA = DREFCO(IA) *AMP3 
AMP2=1.0 

PTIME1=0.0 
PTIME6=GTIMF (TA) 
GSP1=0.0 

GSP6=GSP 


DO 180 IB=1,IA 
PTIME1=PTIME1+PTIME6 
PTIMEG=PTIME1¢PTIME 
GSP1=GSP1+GSP6 
GSP4=GSP1+GSP 
AMPA=AMPA*AMP2*U REFCO(IB) 
AMP1=1.0 


‘PTIME2=0.0 


PTIMES=TIME (IB) 
GSP2=0.0 
GSP5=THICK (IB) *VEL (IB) /VEL (1) 


DO 150 Ic=IB, LAY 
PTIME2=PTIME2+PTIMES 
GSP2=GSP2+GSP5 

GG=0.. 

DO 220 I=1,IC 
GG=THICK (1) *VEL (I) /VEL(1) 
CONTINUE 
GSPF=GG+GSP2+GSP4 
AMP1=DREFCO (IC) *AMP1 
IF(1IG .EQ- 0) GO TO 107 
AFOM=AM PA*AMP1/GSPF 


GO TO 108 


AFOM=AMPA*AMP1 
TFOM=GTIME (IC) +PTIME2+PTIMES 
IND=IND+1 

B=CMPLX (AFOM, TFOM) 

IF (ISUPFO.EQ.0) GO TO 111 
WRITE(6,101) IND 

IBM1=IB-1 
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IF (ISUPFO.EQ.0) GO TO 112 
WRITE(6,102) IA, IBM1,IC 
WRITE(6,103) AFOM 
WRITE(6,104) TFOM 
112 CALL RDOFF(B,P50,IAFG) 

ICP1=Ic+1 
IF(ICP1.GT.LAY) GO TO 150 
AMP1=DTRACO(IC) *UTRACO(ICP1) *AMP1/DREFCO(IC) 
PTIMES=TIME(ICP1) 
GSP5=THICK (ICP1) *VEL(ICP1) /VEL (1) 

150 CONTINUE 


AMP2=1/UTRACO (1B) 

PTIME6=-TIME(IB) 

GSP6=-THICK (IB) *VEL (IB) /VEL (1) 
180 CONTINUE 


TAP1=IA+1 

IF(IAP1 .GT.LAY) GO TO 200 

AMP 3=DTRACO (IA) *UTRACO(IAP1) *AMP3 
200 CONTINUE 


RETURN 
END 


C *SOM* CALCULATES THE ARRIVAL TIME AND AMPLITUDE 
C OF THE SECCND ORDER MULTIPLES 
SUBROUTINE SOM(IND,LAY, TIME, GTIMF, DREFCO, THICK, DEPTH, IIG, 
CDTRACO,UREFCO, UTRACO, ISUPSO,B, P50, I AFG) 
COMPLEX P50(8192) 
DIMENSION GTIME (131) 
COMPLEX B,CMPLX 
DIMENSION DREFCO (131) 
DIMENSION DTRACO (131) 
DIMENSION UREFCO (131) 
DIMENSION UTRACO (131) 
DIMENSION TIME(¢131) 
DIMENSION DEPTH (131) 
DIMENSION THICK(131) 
100 FORMAT(//,T28,'SECOND ORDER MULTIPLES',//) 
101 FORMAT(//,*GROUP SYMBOL! ,T30,15) 
102. FORMAT (*+",T40, 514) 
103. FORMAT('AMP =",£15.5) 
104 FORMAT('TIME =*,F10.5) 
105 FORMAT (F10.5) 
106 FORMAT('GSPF =',E15.4) 
IF(ISUPSO-.EQ.0) GO TO 110 
WRITE (6, 100) 
110 AMP10=1.0 


po..200 IA=1,LAY 
PTIME=GTIME(IA) 
CHIE 
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DO 210 I=1,IA 
GSP=GSP+ THICK (I) *VEL (I) /VEL (1) 
CONTINUE 
AMPA=DREFCO (IA) *AMP10 
AMP9=1.0 

AMP17=AMPA 

PTIME1=0.0 

PTIME6=GTIME(TA) 

GSP1=0.0 

GSP6=GSP 


DO 180 IB=1,IA 
PTIME1=PTIME14+PTIME6 
PTIMEG=PTIME1+PTIME 
GSP1=GSP 1+GSP6 
GSP4=GSP1#GSP 
AMPA=AMP17*AMP9*UREFCO (IB) 
AMP8=1.0 

AMP15=AMPA 

PTIME2=0.0 

PTIMES=TIME(IB) 

6S B2-0.40 
GSP5=THI CK (IB) *VEL (IB) /VEL (1) 


DO 150 IC=IB, LAY 
PTIME2=PTIME2+PTIMES 
GSP2=GSP2+GSP5 
AMP8=DREFCO (IC) *AMP8 
AMP7=1.0 
AMPA=AMP15*AMP8 
AMP16=AMPA 
PTIME3=0.0 
PTIME7=GTIME (IC) 
GSP7=0.0 

GSP8=0. 

DO 42 Osel= 174L0 
GSP8=GSP8+THICK (1) *VEL( I) /VEL(1) 
CONTINUE 

DOw UG OaeD= Wu, 1c 
PTIME3=PTIME3*+PTIMET 
GSP7=GSP7+GSP8 
AMPA=AMP16*AMP7*UREFCO(ID) 
AMP6=1.0 

PTIME8=0.0 

PTIMEQ=TIME (ID) 

GSP9=0.0 
GSP10=THICK (ID) *VEL(ID) /VEL (1) 
DO 130 IE=ID,LAY 
PTIMES=PTIMES*PTIMEY 
GSP9=GSP9+GSP10 
AMP6=AMP6*DREFCO (IE) 

GPPR=0. 

DO RZ spel slice! © 
GPPR=GPPR¢THICK (I) *VEL(I) /VEL (1) 
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130 


140 
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180 


200 


CONTINUE 

GS PF=GPPR+GSP9+GSP7+GSP2+GSPu4 
IF(IIG .EQ.0) GO TO 107 
AFOM=AMPA*AMP6/GSPF 

GO TO 108 

AFOM=AMPA*AMP6 

TFOM=GTIME (IE) ¢PTIME2+PTIMEU+PTIME3+PTIMES 
IND=IND+1 

B=CMPLX (AFOM, TFOM) 

IF(ISUPSO.EQ.0) GO TO 111 

WRITE (6,101) IND 

IBM1=IB-1 

IDM1=ID-1 

IF (ISUPSO.EQ.0) GO TO 112 

WRITE (6,102) IA,IBM1,1C,IDM1,1E 
WRITE(6,103) AFOM 

WRITE(6,104) TFOM 

CALL RDOFF(B,?50,IAFG) 

IEP 1=IEF+1 

IF(IEP1.GT.LAY) GO TO 130 
AMP6=UTRACO (IEP1) *DTRACO (IE) *AMP6/DREFCO(IE) 
PTIMES=TIME(IEP1) 
GSP10=THICK (IEP1) *VEL(IEP1) /VEL(1) 
CONTINUE 

AMP7=1.0/UTRACO(ID) 
PTIME7=-TIME(ID) 
GSP8=-THICK (ID) *VEL(ID) /VEL(1) 
CONTINUE 

ICP1=1C+1 

IF(ICP1.GT.LAY) GO TO 150 
AMP8=UTRACO (ICP1) *DTRACO(IC) *AMP8/DREFCO(IC) 
PTIMES=TIME(ICP1) 
GSPS=THICK (ICP1) *VEL(ICP1) /VEL (1) 
CONTINUE 

AMP9=1.0/UTRACO (IB) 
PTIME6=- TIME (IB) 
GSP6=-THICK (1B) *VEL (IB) /VEL (1) 
CONTINUE 

TAP1=IA+1 

IF(IAP1 .GT.LAY) GO TO 200 
AMP10=DTRACO (IA) *UTRACO (IA+1) *AMP10 
CONTINUE 

RETURN 

END 


C 'psyt PLOTS THE SYNTHETIC SEISMOGRAM 


SUBROUTINE PSY(C,C201, P60) 


C DIMENSION Cc -SEE MAIN 


DIMENSION C (4602) 


C DIMENSION C201( -SEE MAIN 


DIMENSION C201(4602) 
DIMENSION P6041) 
CALL PLOTS 
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MOVE ORIGIN IN TO ALLOW FOR LABELS ETC. 
CALL PLOT(8.0,8.0,-3) 
FOR THE TIME VALUES FORCE 10 INCHES PER SEC. 
XMIN=0.0 
XDELT=.2 
STORE SCALING PARAMETERS FOR THE TIME AXIS 
C201(4601) =XMIN 
C201 (4602) =XDELT 
FOR THE AMPLITUDE AXIS WILL FORCE Te 7LNCH 
YMIN=0.0 
YDELT=1. | 
STORE SCALING PARAMETERS FOR AMPLITUDE AXIS 
C (4601) =YMIN 
C (4602) =YDELT 
PLOT TIME AXIS 
CALLTAXIS(0.0,-1.0,*TIME-SEC!,-8,15.0,0.0,0.0,.2) 
PLOT AMPLITUDE AXIS 
CALL AXIS(0.0,-1.0,*AMPLITUDE',9,2.0,90.0,-1.,1-) 
PLOT DATA 
CALL PLOT(0.0,0.0, 3) 
DO 1 T=1,3000 
XP= (C201 (I-C201 (4601) ) /C201 (4602) 
TP=C (I) /C( 4602) 
CALL PLOT(XP,TP, 2) 
1 CONTINUE 
PLOT SPIKE SYNTHETIC 
MOVE ORIGIN UP 
CALL PLOT(8.0, 16.0, -3) 
PLOT TIME AXIS 
CALtMAKES(0, 0,-1.0, TINE <SEC*/-8,15.070.0,0.0,. 2) 
PLOT AMPLITUDE AXIS 
CALE SAXIS(0.0,-1.0, 'ANPLITUDE",9,2-°50,90.0,-—1., 1.) 
PLOT DATA FOR SPIKE SYNTHETIC 
CALL PLOT(0.0,0.0, 3) 
DO 2 I=1,3000 
XP= (C201 (1-C201 (4601) ) /C201 (4602) 
TP=P60(I) /YDELT 
CALL PLOT (XP,0.0,2) 
CALL. PLOT(XP,TP, 2) 
CALEse PLOT(XP, 0.042) 
CONTINUE 
CALLPPLOT (280,72. 0, 999) 
RETURN 
END 


‘RDOFF' ROUNDS OFF THE TIME OF ARRIVAL OF A SEISMIC PULSE 
"0 3 DECIMAL DIGITS AND PLACES THE AMPLITUDE VALUE IN 
P5O0(IDIV) WHERE IDIV/1000 IS THE TIME OF ARRIVAL TO 
THE CLOSEST MILLISECOND 

SUBROUTINE RDOFF(B,P50,IAFG) 

COMPLEX B 
DIMENSION P50( -SEE MAIN 

COMPLEX P50(8192) 
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TAFG=IAFG+1 

TIME1=( ATMAG(B)) *1000. 

TIME2=TIME1+.5 

IDIV=TIME2 

P50 (IDIV) =CMPLX (REAL (P50 (IDIV) ) #+REAL(B) ,0.0) 
RETURN 

END 


SUBROUTINE NLOGN(N,X,SIGN) 

C DIMENSION MMM ¢ -SEE MAIN 
DIMENSION MMM (13) 

C DIMENSION X( -SEE MAIN 
DIMENSION X(8192) 
COMPLEX X,WK,HOLD,Q 
LX=2**N 
Dory lad yN 

1 MMM (I) =2** (N-T) 

DO 4 L=1,N 
NBLOCK=2 ** (L=1) 
LBLOCK=LX/NBLOCK 
LBHALF=LBLOCK/2 
K=0 
DO 4 IBLOCK=1, NBLOCK 
FK=K 
FLX=LX 
V=SIGN¥*6.2831853*FK/FLX 
WK=CMPLX (COS (V) ,SIN(V)) 
ISTART=LBLOCK* (I BLOCK-1) 
DO 2 I=1,LBHALF 
J=ISTART+I 
JH=J+LBHALF 
Q=X (JH) *WK 
X (JH) =X (J-Q 
X(J) =X (J) +Q 

2 CONTINUE 
DO 3 I=2,N 
II=1I 
IFq(Ke-LT.MMMq(I)) GO TO 4 

3 K=K-MMM(T) 

4 K=K+MMM (II) 
K=0 
DO 7 J=1,LX 
IF(K.LT.J) GO TO 5 
HOLD=X (J) 
X (J) =X (K+1) 
X (K+ 1) =HOLD 

5 po 6 I=1,N 

II=1I 

IF (Ke-LT.MMM(I)) GO TO 7 

K=K-MMM (qT) 

K=K+MMM (ITI) 

IF(SIGN.LT.0.0) RETURN 

po 8 i{=1,LX 
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8 X(T)=X (1) /FLX 
RETURN 
END 


SUBROUTINE KLAUD GENERATES A NORMALIZED KLAUDER WAVELED 
F1 - THE TERMINAL LOW FREQUENCY 
IZ - THE FREQUENCY BANDWIDTH IN OCTAVES 
T - THE DURATION OF THE SIGNAL IN SECONDS 
THE GENERATION IS DONE USING EQN.7 S.S.C. TRAINING DEPT. 
CATALOG OF KLAUDER WAVELETS 
SUBROUTINE KLAUD (T,F1,Z,VALUE) 
DIMENSION VALUE( ) IN THE FOLLOWING STATEMENT-SEE MAIN 
COMPLEX VALUE(8192), CMPLX 
PI=3.141592 
"VALUE (301) IN THE NEXT STATEMENT REPRESENTS THE AUTO- 
CORRELATION OF THE SWEEP AT ZERO SHIFT AND SHOULD 
THEORETICALY HAVE A VALUE OF'T', WHERE 'T* IS THE 
TIME LENGTH OF THE SWEEP. IF NORMALIZATION IS 
DESIRED SPECIFY "COMPLEX(1.0,0.0)* RATHER THAN 
COMPLEX(T,0.0)IN THE NEXT STATEMENT 
VALUE (301) =CMPLX (1.0, 0.0) 
T=301 
SHIFT=.001 
DO 16 J=1,300 
A=PI*F1*SHIFT* ( (2**Z-1) * (1-SHIFT/T) 
B=PI*F1* SHIFT *((2**Z=1) 
C=PI*F1*SHIFT *( (2**Z) #1) 
VALUE (I+J) =CMPLX (1. 0*SIN (A) *COS (C) /B, 0.0) 
VALUE (I-J) =VALUE(I+J) 
SHIFT=SHIFT+.001 
16 CONTINUE 
DO 18 I=602,8192 
VALUE (I) =CMPLX(0.0,0.0) 
18 CONTINUE 
RETURN 
END 
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C SUBROUTINE WAVLT1 STORES IMP. LOG WAVELET 
C IN VALUE 
SUBROUTINE WAVLT1 (VALUE) 
COMPLEX VALUE (8192) ,CMPLX 
DIMENSION PALL(97) 
10 FORMAT (5F10.5) 
READ (13,10) (PALL(1I),I=1,97) 
pO 15 I=1,97 
VALUE (I) =CMPLX (PALL(1) , 0.0) 
15 CONTINUE 
DO 20 I=98,8192 
VALUE (I) =CMPLX (0.0,0-0) 
20 CONTINUE 
RETURN 
END 
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C.2 The Permuted Partition Approach 


C "FACT" IS CALLED BY *NODYAN* AND COMPUTES NJ FACTORIAL 
SUBROUTINE FACT (AJ,AJF) 
IF (AJ .EQ. 0.0) GO TO 14 
IF (AJ .EQ. 1.0) Go TO 14 
AJM2=AJ-2.0 
NIM2=AIM2 : 
IF(AJM2 .FQ. 0.0) GO TO 15 
AJF=AJ 


DO 13 I=1,NJIM2 
BJ=AJ-I 
AJF=AJF*BJ 

13 CONTINUE 


GO TO 16 
14 AJF=1.0 
GO TO 16 
15 AJF=AJ 
16 RETURN 
END 


C "NODYAN CALCULATES THE NUMBER OF DYNAMIC ANALOGS 
C ASSOCIATED WITH A PARTICULAR RAY CODE 
C REFERENCE EQN. 19 HRON 1972 

SUBROUTINE NODYAN(MAC,ICN,ICM, AF) 

DIMENSION ICN(17) 

DIMENSION ICM (16) 

IC=MAC-1 

AF=1.0 


DO 18 T=1,IC 
A15=ICN (I) 
CALL FACT(A15,AJF) 
A20=ICM (I) 
CALL FACT(A20,AJG) 
CALL FACT((A15-A20) ,AJH) 
AZ1=AJF/ (AJG*AJH) 
IP1=I+1 
A30=ICN (IP1) 
CALL FACT((A30-1.0),AJT) 
A60=A30-A15+A20 
CALL FACT(A60,AJJ) 
ASO=A15-A20-1.0 
CALL FACT((A50),AJK) 
AZ2=AJI/ (AJJ* AIK) 
ANC=AZ1*AZ2 
AF=AF*ANC 

18 CONTINUE 
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RETURN 
END 


"RDOFF* ROUNDS OFF THE TIME OF ARRIVAL OF A SEISMIC PULSE 
TO 3 DECIMAL DIGITS AND PLACES THE AMPLITUDE VALUE IN 
PSO(IDIV) WHERE IDIV/1000 IS THE TIME OF ARRIVAL TO 
THE CLOSEST MILLISECOND 

SUBROUTINE RDOFF(B,P50) 

COMPLEX B 

COMPLEX P50(8192) 

TIME1= (AIMAG (B) ) *1000. 

TIME2=TIME1+.5 

IDIV=TIME2 

P50 (IDIV) =CMPLX (REAL (P50 (IDIV) ) #REAL (B) ,0.0) 

RETURN 

END 


SUBROUTINE NLOGN(N,X,SIGN) 
DIMENSION MMM (13) 
DIMENSION X (8192) 

COMPLEX X,WK,HOLD, 0 
LX=2**N 


pO 1 I=1,N 
1 MMM «(Ip =2%* (N-T) 


DO 4 L=1,N 
NBLOCK=2 ** (L-1) 
LBLOCK=LX/NBLOCK 
LBHALF=LBLOCK/2 
K=0 


DO 4 IBLOCK=1,NBLOCK 
FK=K 

FLX=LX 
V=SIGN*6. 28 31853*FK/FLX 
WK=CMPLX (COS (V) ,SIN(V)) 
ISTART=LBLOCK* (IBLOCK -1) 


DO 2 I=1,LBHALF 
J=ISTART+I 
JH=J+ LBHALF 
Q=X (JH) *WK 
X (JH) =X (J-Q 
X (J) =X (J) + 

2 CONTINUE 
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79 
3 K=K-MMM (I) 
4 K=K+MMM(II) 
K=0 


DOF amdesle, LX 
IF(K.LT.J) GO TO 5 
HOLD=X (J) 

X (J) =X (K+1) 

X (K+1) =HOLD 


5 enrages —1, N 

II=1I 

IF(K.LT.MMM(I)) GO TO 7 
6 K=K-MMM(T) 


7 K=K+#MMM (IT) 
IF (SIGN.LT.0.0) RETURN 


DO 8 I=1,LX 
8 X(T) =X(I)/PFLX 


RETURN 
END 


C "TAMP* CALCULATES THE TIME OF ARRIVAL AND AMPLITUDE OF A 
C PARTICULAR PHASE GIVEN ONLY THE RAY CODE 
SUBROUTINE TAMP (ICP1,LAY,ICM,ICN, DREFCO, DTRACO,IIG, THICK, 
CUREFCO, UTRACO, TIME, AF, B,P50, VEL) 
COMPLEX B,CMPLX 
COMPLEX P50(8192) 


50 
51 


DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 


FORMAT (*AMP =",E15.4) 
FORMAT(* TIME =',F10.5) 


ICM (ICP1) = 
AMP=1.0 


I3(16) 
ICN(17) 
ICM(16) 
DREFCO (16) 
DTRACO (16) 
UREFCO (16) 
UTR ACO (16) 
TIME (16) 
THICK (16) 
VEL(15) 


ICN (ICP1) 


t=1 

GSPF=0.0 
TIM=0.0 
T3¢1) =1 


DO 92) J=27UAY 
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100 


101 


I3 (J) =0 
CONTINUE 


TEQICOCLT scl... 0) GOTO 54 
AMP=AMP*DTRACO(T) 

I=I+1 

I3 (I) =13 (1) +1 


GO TO 53 


IND50=1 

AMP=AMP*DREFCO(I 

IM1=I-1 | 

IF(IM1 .EQ. 0) GO TO 100 


DO 56 IB=1,1M1 

IMIB=I-IB 

IZ=IB 

IF(ICNqIMIB) .EQ. I3(IMIB)) GO TO 59 
TIB1=I-IB+tl1 

AMP =AMP*UTRACO(IIB1) 

CONTINUE 


AMP=AMP*UREFCO(IMIB) 

IF(IND5O .EQ. ICM(I)) GO TO 60 
IND5O0 =IND50+1 

IM1=I-1 


DOW oT LO, LMI 
AMP=AMP*DTRACO (IB) 
CONTINUE 


DO 58 J=1,1 
I3 (J) =13 (J) +1 
CONTINUE 


GO TO 55 


IF(IND50 .EQ. ICMqI)) GO TO 101 
AMP=AMP*UREFCO (1) 

13 (I) =13(I) +1 

IND5O=INDSO+1 


GO TO 55 


IF(I .EQ.ICP1) GO TO 150 
AMP=AMP*UREFCO (1) 

13 (1) =13 (1) +1 
AMP=AMP*DTRACO(1) 

T=I+1 

13(1) =13 q1) +1 


GO TO 53 
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59 IF(INDSO .EQ. ICM(I)) GO To 140 
IIB1=I1-1I2Z+1 
AMP=AMP*URFFCO(IIB1) 
IND50=IND50+1 
I3 (1) =13 (1) +1 


GO TO 55 

60 "*D0T104, 34=13 7 
AMP=AMP*DTRACO(I1) 

104 CONTINUE 
IP1=I+1 
DO 105 J=1,IP1 
13 (J) =13 (J) +1 

105 CONTINUE 
I=I+1 
GO TO 53 

140 IF(I .EQ. ICP1) GO TO 150 
TIB1=I-IZ+1 
AMP=AMP*UREFCO(IIB1) 
DO 141 J=IIB1,1 
AMP=AMP*DTRACO (J) 

141 CONTINUE 
IP 1=1+1 
DO 142 J=IIB1,1P1 
3 (J) =13 (J) +1 

142 CONTINUE 
T=I+1 
GOe TOS 53 

150 DO 151 IT=1,ICP1 
AMP=AMP*UTRACO (I) 

151 CONTINUE 
IF (IIG .EQ. 0) GO TO 154 
DO 153 I=1,ICP1 
GSPF=GSPF+t2* ICN (I) *THICK (I) *VEL (I) /VEL(1) 

153 CONTINUE 
AMP=AMP*AF/GSPF 


GO TO 155 
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C SPECIFY THE TIME COTOFF IN THE NEXT STATEMENT 


Cc 
Cc 


154. 
155 


152 


AMP=AMP*AF 


BUI Zi =1. TCD 4 
TIM=TIM+ 2* ICN (I) *T IME (I) 
CONTINUE 


TE(TIM CGT. oe. )8 GO TC 257 
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MOVE ORIGIN IN TO ALLOW FOR LABELS ETC. 


FOR THE TIME VALUES FORCE 5 INCHES PER SEC. 


STOKE SCALING PARAMETERS FOR THE TIME AXIS 


FOR THE AMPLITUDE AXIS WILL FORCE 1./INCH 


STORE SCALING PARAMETERS FOR AMPLITUDE AXIS 


PLOT 
PLOT 


PLOT 


3 
PLOT 


WRITE (6,50) AMP 

WRITE(6,51) TIM 

B=CMPLX (AMP, TIM) 

CALL RDOFF(B, P50) 
RETURN 

END 


PLOTS THE SYNTHETIC SEISMOGRAM 
SUBROUTINE PSY(C,C201, P60) 
DIMENSION C(4602) 

DIMENSION P60 (4602) 

DIMENSION €C201(4602) 

CALL PLOTS 


CALL PLOT(8.0,8.0,~3) 


XMIN=0.0 
XDELT=.2 


C201( 4601) =XMIN 
C201(4602) =XDELT 


YMIN=0.0 
YDELT=.02 


C (4601) =YMIN 

C (4602) =YDELT 
P60(4601) =YMIN 
P60 (4602) =YDELT 
TIME AXIS 


CALL AXIS (00 0p leGe (LOS pO se Vis Uisle gine) 


AMPLITUDE AXIS 


CALL AXIS(0.0,-1.0,°AMPLITUDE',9,2.0,90.0, 


DATA 
CALL PLOT(0.0,0.0, 3) 

po 1 T=1,200 

XP= (C201 (I-C201 (4601) ) /C201 (4602) 
TP=C (1) /C( 4602) 

CALL PLOT(XP,TP, 2) 

CONTINUE 

pO 3 1=201,3000 

XP=C201(1) 

TP=C(I) /-0001 

CALL PLOT(XP,TP, 2) 

CONTINUE 

SPIKE SYNTHETIC 
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MOVE ORIGIN UP 
CALL PLOT(8.0,16.0,-3) 
PLOT TIME AXIS 


CALL AXTS(0.0,-1.0,"TIME-SEC',-8,15. TAs Oe Oy teh 
PLOT AMPLITUDE AXIS 


CALL AXTS(0.0,-1.0, AMPLITUDE" ,9,2.0,90.0,-12, 1) 
PLOT DATA FOR SPIKE SYNTHETIC 
CALL PLOT(0.0,0.0, 3) 
PQ y2)1=1,200 
XP= (C201 (I-C201( 4601) ) C201 (4602) 
TP=P60 (IT) /P60 (4602) 
CALL PLOT (XP,0.0,2) 
CALL PLOT(XP,TP, 2) 
CALL PLOT(XP,0.0,2) 
CONTINUE 
DO 4 T=201,3000 
XP=C201 (1) 
TP=C(I) 7.0001 
GARG PEOTEXP,0.,2) 
CALL PLOT(XP,TP, 2) 
CALL PLOT(XP,0.,2) 
4 CONTINUE 
CALL PLOT(2.0,2. 0,999) 
RETURN 
END 


SUBROUTINE KLAUD GENERATES A NORMALIZED KLAUDER WAVELET 
F1 - THE TERMINAL LOW FREQUENCY 
IZ - THE FREQUENCY BANDWIDTH IN OCTAVES 
T - THE DURATION OF THE SIGNAL IN SECONDS 
THE GENERATION IS DONE USING EQN.7 S.S.C. TRAINING DEPT. 
CATALOG OF KLAUDER WAVELETS 
SUBROUTINE KLAUD (TyF1,2, VALUE) 
DIMENSION VALUE( ) IN THE FOLLOWING STATEMENT-SEE MAIN 
COMPLEX VALUE(8192), CMPLX 
PT=3.141592 


'VALUE(301) IN THE NEXT STATEMENT REPRESENTS THE AUTO- 
CORRELATION OF THE SWEEP AT ZERO SHIFT AND SHOULD 
THEORETICALY HAVE A VALUE OF'T', WHERE 'T* IS THE 
TIME LENGTH OF THE SWEEP. IF NORMALIZATION IS 
DESIRED SPECIFY "COMPLEX(1.0,0.0)* RATHER THAN 
COMPLEX (T,0.0)IN THE NEXT STATEMENT : 

VALUE (30 1) =CMPLX (1.0,0.0) 

T=301 

SHIFT=.001 


DO 16 J=1,300 
A=PI*FI*SHIFT* ((2**Z- 1) * (1-SHIFT/T) 
B=PI*F1* SHIFT *( (2**Z-1) 
C=PI*F1*SHIFT *( (2**Z) +1) 

VALUE (I+d) =CMPLX (1.0*SIN (A) *COS (C) /B, 0.0) 
VALUE (I-J) =VALUE(I+J) 
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C SUBROUTINE WAVLT1 STORES IMP. 


SHIFT=SHIFT+.001 


DO 18 I=602,8192 


VALUE (I) =CMPLX (0.0,0.0) 


16 CONTINUE 
18 CONTINUE 
RETURN 
END 
LUE 


C IN VA 


10 


15 


20 


SUBROUTINE WAVLT1(VALUE) 
COMPLEX VALUE (8192) ,CMPLX 


DIMENSION 


PALL (97) 


PORMAT (SF10.5) 
(PALL(I) , T=1, 97) 


READ(7, 10) 


DO 15 I=1, 


97 


LOG WAVELET 


VALUE (I) =CMPLX (PALL(I) , 0.0) 


CONTINUE 


DO 20 T=98,8192 


VALUE (I) =CMPLX (0.0,0.0) 


CONTINUE 


RETURN 
END 
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SUBROUTINE PERMUT(MAC,ICN, ICM, MAXLAY,IND,LAY,DREFCO, 
CDTRACO ,IIG, THICK, UREFCO, UTRACO,TIME,AF, IT AFG,B,P50,VEL) 


DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSTON 
DIMENSION 
DIMENSION 
DIMENSION 


N (16,16) 
ICN (17) 
ICM (16) 
IBEN (16) 
DREFCO (15) 
DTRACO (16) 
THICK (16) 
UREFCO (16) 
UTRACO (16) 
TIME (16) 
VEL(15) 


COMPLEX P50(8192) 


DIMENSION 
COMPLEX B 


J (16) 


ITT=16-MAC+1 


pO 10 J9=1,16 


DO 9 K=ITT, 16 
N(J9,K) =ICN(K-ITT+1) 
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9 CONTINUE 
10 CONTINUE 


CALL DICM(ICN,ICM,MAXLAY,MAC,IND,LAY,DREFCO, 
CDTRACO ,IIG, THICK, UREFCO, UTRACO, TIME, AF, 1AFG,B,P50,VEL) 


Om m= 1) 16 
Tet 
11 CONTINUE 


GO TO 1000 


800 IF(J(J2).LE. J3) GO TO 806 
IM1=J (J 2-1 


DO 802 I=J3,JM1 
IF(N(J2,J3(J32)).EQ.N(J2,1)) GO TO 1004 
802 CONTINUE 


806 IF (N(J2,32).FQ.N(J2,J(J2))) GO TO 1004 


pO 801. I=1,MAC 

ITMAC=16-MAC+4I 

N(J3,ITMAC) =N(J2,ITMAC) 
801 CONTINUE 


N(J3,52) =N (J2,J5 (J2)) 
N(J3,J (J2) ) =N(J2,J2) 


pO 805 I=1,MAC 
IBEN (I) =N(J3, 16-MAC#TI) 
805 CONTINUE 


CALL DICM{IBEN,ICM,MAXLAY,MAC,IND,LAY,DREFCO, 
CDTRACO ,I1G, THICK, UREFCO, UTRACO, TIME, AF, IAFG,B,P50,VEL) 
ISEV=16-J1+3 


DO 804 L=ISEV, 16 


DO 803 Kk=1,16 
N(L,K)=N(J3,K) 
803 CONTINUE 


8Qu CONTINUE 


1000 IF (N(15, 15) .EQ.N(15,16)) GO TO 1001 
N(16, 15) =N (15,16) 
N (16, 16) =N (15,15) 


po 904 T=1,MAC 
IBEN (I) =N (16, 16-MACtT) 
904 CONTINUE 
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CALL DICM(IBEN,ICM,MAXLAY,MAC,IND,LAY,DREFCO, 
CDTRACO ,ITIG, THICK, UREFCO,UTRACO, TIME, AF, LAFG,B, P50,VEL) 


GO TO 1001 
1004 IFR=J1 

GO TO 1003 
1001 IFR=3 


1003 DO 1010 J1=IFR, 16 
IF(MAC.LT.J1) GO TO 1020 
J2=16-J1+1 
J (52) =J3 (32). +1 
J3= 16-3142 


DO 1002 J4=J3,16 
J(J4) =J4 
1002 CONTINUE 


IF(J(J2) -LE.16) GO TO 800 
1010 CONTINUE 


1020 RETURN 
END 


C DICH DETERMINES ,THE VECTOR ICM 
SUBROUTINE DICM(ICN, ICM, MAXLAY,MAC,IND,LAY,DREFCO, DTRACO, 
CLIG, THICK,UREFCO, UTRACO, TIME, AF,IAFG,B,P50,VEL) 


DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 


MMN (16) 
MMX (16) 

M (16) 

ICN (17) 
ICM (16) 
DREFCO (16) 
DTRACO (16) 
THICK (16) 
UTRACO (16) 
UREFCO (16) 
TIME (16) 
VEL (15) 


COMPLEX P50 (8192) 


COMPLEX B 


23 FORMAT(//, ‘GROUP NUMBER',T25,15) 
24 FORMAT(*ICN =',725,1015) 
25 FORMAT (*ICM =",725,1015) 
26 FORMAT(*N =',F8. 3) 
IC=MAC-1 


DO 260 I=1,1C 
MMX (I) =ICN(I-1 
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MMN (I) =MAXO (ICN (I- ICN (I+#1), 0) 
260 CONTINUE 


264 M(1) =MMN (1). 
I25=MAXLAY-2 
309 I20=1 


30 (DOL 35061 15=120, 725 
IF (IC. EQ.115) GO TO 399 
T30=115+1 
M (130) =MMN (130) 

390 CON TUNILE 


BERS) A O52 at 


400 DO 450 I=1,IC 
ICM (I) =M (1) 
450 CONTINUE 


IND=IND+1 
WRITE(6,23) IND 
WRITE (6,24) (ICN(I),I=1,ICP1) 
WRITE (6,25) (ICM qI),1I=1,1C) 
CALL NODYAN (MAC, ICN,ICM, AF) 
WRITE(6,26) AF 
IJAF=IFIX (AF) 
ITAFG=IAFGtIJAF 
CALL TAMP(MAC,LAY,ICM,ICN, DREFCO, DTRACO,IIG, THICK, UREFCO, 
CUTRACO,TIME,AF,3B,P50,VEL) 
331 I200=MAXLAY-2 
MLM1=MAXLAY=-1 


DO 520 1100=1,1200 

I201=MLM1-1100 

CL BOs 20 1) G0 oO. og 
520 CONTINUE 


M(MLM1)=M(MLM1) +1 
IF(M(MLM1) «LE. MMX(MLM1)) GO TO 400 
559 IF (1201.EQ.1) GO T01011 
560 M1201) =M (1201) +1 
IF(M(I201) .LE. MMX(I201)) GO TO 1020 
1201=1201-1 


GO TO 559 
1011 M™q1)=M(1) +1 

IF(M(1) .LE. MMX(1)) GO TO 309 

IF (M(1) .GT. MMX(1) ) GO TO 650 
LOZOmar eC E20 

GO TO 310 


650 RETURN 
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END 
MAIN FOLLOWS 


500 FORMAT (//,T1C0,*DOWNWARD REFLECTION COEFFICIENTS',//) 
501 FORMAT (//,T10,"DOWNWARD TRANSMISSION CORRPICIENTS* 7/7) 


502 FORMAT (//,T10,*UPWARD REFLECTION COEFFICIENTS! ,//) 


503 FORMAT(//,T10,*UPWARD TRANSMISSION COEFFICIENTS',//) 
14 FORMAT (F10. 2) 
9  FORMAT(8F10.5) 
10 FORMAT (8F10.5) 
11. +FORMAT(8F10.5) 
12 FORMAT(8F10.5) 
DIMENSION VEL(NO. OF VELOCITY VALUES AVAILABLE) 
DIMENSION VEL (15) 
DIMENSION DEPTH(NO. OF VELOCITY VALUES - 1) 
DIMENSION DEPTH (14) 
LAY = NO. OF VELOCITY VALUES - 1 
LAY=14 
THICK (I) CONTAINS THE THICKNESS OF THE (1)TH LAYER 
THE VALUES THICK(I) ARE CALCULATED INTERNALLY GIVEN THE 
DEPTH VALUES 
DIMENSION THICK (NO. OF VELOCITY VALUES - 1) 
DIMENSION THICK(14) 
DREFCO REFERS TO THE DOWNWARD REFLECTION COEFFICIENT 
DIMENSION DREFCO(NO. OF VELOCITY VALUES - 1) 
DIMENSION DREFCO(14) 
DTRACO REFERS TO THE DOWNWARD TRANSMISSION COEFFICIENT 
DIMENSION DTRACO(NO. OF VELOCITY VALUES - 1) 
DIMENSION DTRACO(14) 
UREFCO REFERS TO THE UPWARD REFLECTION COEFFICIENT 
DIMENSION UREFCO(NO. OF VELOCITY VALUES - 1) 
DIMENSION UREFCO (14) 
UTRACO REFERS TO THE UPWARD TRANSMISSION COEFFICIENT 
DIMENSION UTRACO(NO. OF VELOCITY VALUES - 1) 
DIMENSION UTRACO (14) 
TIME(I) CONTAINS THE ONE WAY VERTICAL TRAVEL TIME THROUGH 
LAYER I 
DIMENSION TIME(NO. OF VELOCITY VALUES - 1) 
DIMENSION TIME(14) 
VALUE CONTAINS THE KLAUDER (RIKER) WAVELET PLUS TRAILING 
ZEROS. THE DIMENSION OF VALUE MUST BE GREATER THAN OR 
EQUAL TO THE LENGTH OF THE CONVOLUTION OF THE WAVELET 
WITH THE SPIKF SYNTHETIC, WITH THR ADDED CONDITION THAT 
THE DIMENSION MUST EQUAL 2**N WHERE N IS SOME POSITIVE 
INTERGER. 2e*5=32 2**6=64 2**7=128 2**8=250 2**9=512 
2**10=1024 2**11=2048 2**12=4096 2**13=8192 2**14=16384 
COMPLEX VALUE(8192) 
P50 CONTAINS THE SPIKE SYNTHETIC (OUTPUT FROM RDOFF) 
DIMENSION PSO(SAME DIMENSION AS VALUE) 
COMPLEX P50(8192) 
DIMENSION BNUIS(2 * DIMENSION OF VALUE) 
DIMENSION BNUIS (16384) 
EQUIVALENCE (VALUE, BNUIS) 
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C CONTAINS THE Y OR AMPLITUDE VALUES USED IN PLOTING THE 
SYNTHETIC SEISMOGRAM 
DIMENSION C(LENGTH OF THE CONVOLUTION OF THE SPIKE SYNTHETIC 


WITH THE WAVELET + 2) THE SCALING PARAMETERS ARE STORED IN 
THE LAST TWO POSITIONS 


DIMENSION C(4602) 
C201 CONTAINS THE X OR TIME VALUES USED IN PLOTTING THE SS 
DIMENSION C201(SAME AS DIMENSION OF C) 

DIMENSION C201(4602) 
IN NLOGN THE DATA LENGTH MUST BE EQUAL TO 2**N WHERE N IS 
SOME POSITIVE INTERGER. DIMENSION MMM(LARGEST VALUE OF N TO 
BE PROCESSED) 

DIMENSION MMM (13) 
DIMENSION 13 (LAY) 

DIMENSION 13(14) 

DIMENSION IBEN(14) 

DIMENSION P60(4600) 
IF THE GEOMETRICAL SPREADING FACTOR IN AMPLITUDE IS NOT 
DESIRED SPECIFY IIG=O. ANY OTHER VALUE FOR IIG WILL 
INVOKE THE SPREADING CALCULATION. BE CERTAIN TO 
CHANGE THE VERTICAL SCALING PARAMETERS IN THE 
PLOTTING ROUTINE WHEN RUNNING THIS PROGRAM 
WITHOUT GEOMETRICAL SPREADING. 

TIG=1 

UTRACO (1)=1.0 

UREFCO(1)=1.0 
INPUT VELOCITY AND DEPTH VALUES 

READ(8,14) (DEPTH(I) ,I=1,LAY) 

JISO=LAY+1 

READ(10,14) (VEL (I),I=1,JJ50) 
INSERT DO LOOP PARAMETER ...1=1,DIMENSION OF VALUE 


po 13 I=1,8192 
P50 (I) =CMPLX(0.0,0.0) 
13 CONTINUE 


THICK (1) =DEPTH(1) 
LM1= LAY-1 


DOnce =o, LAY 

IM1=I-1 

THICK (1) =DEPTH (I-DEPTH(IM1) 
8 CONTINUE 


NB epeS osc Ea hed BIE 
TIME (I) =THICK (I) /VEL (I) 
3. CONTINUE 


pO 4 I=1,LAY 
IPp1=I+1 
DREFCO (I) = (VEL (I-VEL (IP1)) / (VEL(I) +VEL (IP1)) 


4 CONTINUE 


WRITE(6, 500) 
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WRITE (6,9) (DRE¥CO(I) ,I=1,LAY) 


DO 5 I=1,LAY 
DTRACO (I) =1+#DREFCO (I) 
5 CONTINUE 


WRITE(6,501) 
WRITE(6,10) (DTRACO(I) ,I=1, LAY) 


DO 6 I=2,LAY 

IM1=I-1 

UREFCO (1)= (VEL (I-VEL(IM1)) /(VEL (IM1) #+VEL(I)) 
6 CONTINUE 


WRITE (6,502) 
WRITE(6,11) (UREFCO(I),I=1,LAY) 


DO GACE=2), LAY 
UTRACO (I) =1. +UREFCO(I) 
7 CONTINUE 


WRITE (6, 503) 

WRITE(6,12) (UTRACO(I) ,1=1,LAY) 
AFTER EACH RAY CODE IS GENERATED THE AMPLITUDE AND TIME OP 
ARRIVAL ARE CALCULATED AND STORED IN B. 

COMPLEX B 
DIMENSION MMN(MIN(IHMS,LAY)) 

DIMENSION MMN(14) 

DIMENSION MMX (SAME AS MMN) 

DIMENSION MMX(14) 

DIMENSION ICN(SAME AS MMN + 1) 

DIMENSION ICN(15) 

DIMENSION ICM(SAME AS MMN) 

DIMENSION ICM(14) 
DIMENSION M(SAME AS MMN) 

DIMENSION M(14) 

MAXLAY=LAY 

IND=0 

TAFG=0 
INPUT IHMS = HALF THE MAX. NO. OF SEGMENTS ALLOWED 

THMS=14 

21 FORMAT (//,*GROUP NUMBER',T25,15) 
22 FORMAT(//,T31,°FINISHED',//) 

23 FORMAT(//,*GROUP NUMBER',T25,15) 
24 FORMAT(*ICN =*,T25,10I5) 

25 FORMAT(*ICM =",1725,1015) 

26 FORMAT(*N =',F10.5) 

27. FORMAT(*N = 1°) 

50 FORMAT(*AMP =',F10.5) 

51 FORMAT('*TIME =',F10.5) 

28 FORMAT(*ICN =*,225,15) 

504 FORMAT(/,*TOTAL NUMBER OF PHASES',1T30,110) 
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C PARTITION 


3000 FORMAT('*MADE IT*) 
4000 FORMAT('PARTITION =",125,1015) 


DO 100 N=1,IHMS 

MQ=1 
IN THE NEXT STATEMENT "IF(N.GT.X).. 
TO OBTAIN THE UNRESTRICTED PARTITIONS SET 
X > HMS. FOR SEVERITY LEVEL 1, X=5...FOR 
SEVERITY LEVEL 2, X=4...FOR SEVERITY LEVEL 3 
X=3...FOR SEVERITY LEVEL 4, X=2. 

IF(N.GT.15) GO TO 85 

GO To 71 
C IN THE NEXT STATEMENT "MQ=N-Y"...FOR SEVERITY 
C LEVEL 1 SET Y=4...SEVERITY LEVEL 2, Y=3... 
C SEVERITY LEVEL 3, Y=2...SEVERITY LEVEL 4, Y=1 
85 MQ=N-2 

MAC=MQ 

GO TO 70 
71 IND=IND+1 

MAC=MQ 

ICN (1) =N 

WRITE (6,23) IND 

WRITE (6,28) ICN(1) 

WRITE(6,27) 

IAFG=IAFG+1 

CALL TAMP(MAC,LAY, ICM, ICN, DREFCO,DTRACO,IIG, THICK, 

CU REFCO, UTRACO, TIME, 1.0,B, P50, VEL) 


(@) (ie) (e) (2) 


Oren@ 


GO TO 80 
70 MACM1=MAC-1 


DO 72 I=1,MACM1 
ICN (I) =1 
72 CONTINUE 


74 TSUM=ICN (1) 
IF (MACM1.EQ.1)GO TO 78 


DO 76 I=2,MACM1 
ISUM=ISUM4ICN ¢I) 
76 CONTINUE 


78 ICN (MAC) =N-ISUM 
CALL PERMUT (MAC, ICN, ICM, MAXLAY,IND,LAY, DREPCO,DTRACO, 


CIIG, THICK ,UREFCO, UTRACO,TIME,AF, TAFG,B,P50,VEL) 
80 IT=MaAc-1 
82 IF (IT.EQ.0) GO TO 88 
IF((ICN(MAC~ICN(IT)).GT.1)GO TO 84 


IT=IT-1 


GO TO 82 
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84 ICRAP=ICN(IT) 


DO 86 T=IT,MACM1 
ICN (1) =ICRAP+1 
86 CONTINUE 


GO TO 74 


88 MAC=MACt#1 
IF(MAC.LE.N)GO TO 70 
100 CONTINUE 


2000 WRITE(6,504) IAFG 
WRITE (6, 22) 
C GENERATE THE TIME AXIS FOR PLOTTING 
C201(1)=0. 
C INSERT DO LOOP PARAMATER. I=2,LENGTH OF TRACE 
C IN SECONDS X 1000 


DO 2001 I=1,4599 
_C€201(1+1) =1/1000. 
2001 CONTINUE 


WRITE(7, 2004) (P50 (1) ,1=1,8192) 
FORMAT (2F10. 5) 

SPECIFY THE KLAUDER WAVELET PARAMETERS 
CALL KLAUD(7.,6.,4-0, VALUE) 
CALL WAVLT1(VALUE) 


(ei (a) (2) © 


DO 2004 TI=1,4602 
POO (I) =REAL(P50(7T) ) 
2004 CONTINUE 


WRITE(11) P60 
C IN THE NEXT TWO STATEMENTS THE FIRST PARAMETER OF 
C NLOGN MUST EQUAL THE DIMENSION OF MMM( ) 

CALL NLOGN(13,VALUE,-1.0) 

CALL NLOGN(13,P50,-1.0) 


DO 2002 I=1,8192 
VALUE (1) =VALUE (I) * P50 (I) 
2002 CONTINUE 


C INSERT 1ST NLOGN PARAMETER 

CALL NLOGN(13,VALUE,1.0) 
C INSERT DO LOOP PARAMETER. I=1,LENGTH OF TRACE 
C IN SECONDS X 1000 


DO 2003 I=1, 4600 
C(I) =BNUIS(2*I-1) 
2003 CONTINUE 
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